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Abstract 

iTiis  thesis  advances  boundary  element  techniques  to  model  thermal  oxidation  of  silicon  in 
two  dimensions.  At  temperatures  encountered  in  thermal  oxidation,  silicon  dioxide  flows 
viscoelastically.  A  reduced-dimension,  generalized  boimdaiy  element  method  for  modeling 
such  a  problem  has  been  developed.  With  a  Laplace  transform  technique,  a  viscoelastic 
kernel  function  in  derived  from  Kelvin’s  solution,  which  is  the  fundamental  solution  to 
linear  elasticity.  Constant-velocity  loading  is  chosen  to  operate  with  a  wide  range  of  stress 
relaxation  times.  This  scheme  is  capable  of  replacing  boundary  element  methods 
developed  for  slow  viscous  flow  and  elastic  deformation.  The  oxidant  diffusion  problem  is 
solved  using  a  standard  potential  method  for  Laplace  problems.  Generated  by  oxide 
growth,  stress  affects  both  oxidant  diffusion  and  oxide  flow.  In  particular,  it  changes  the 
diffusivity  of  oxidants  and  viscosity  of  oxide,  rendering  the  diffusion  and  flow  problems 
nonhomogeneous.  Domain  solutions  are  needed  for  both  processes.  A  unified  initial 
stress/built-in  field  formulation  has  been  developed  to  account  for  the  nonlinear  effects, 
using  interior  cells  that  are  placed  where  stress  is  significant.\The  interior  solutions  are 
realized  with  an  interfacial  source  method,  whereby  an  area  i^egral  for  a  cell  is 
transformed  to  a  line  integral  on  the  perimeter  of  the  cell.  It  h^  been  found  that  kernel 
functions  based  on  Kelvin’s  solution  are  deficient  in  modeling  incompressible  materials 
with  a  “hole”.  A  correction  method  that  uses  a  source  term  ^aced  at  the  center  of  the  hole 
has  been  implemented  to  overcome  the  numerical  problenu 
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ABSTRACT 


This  thesis  advances  boundary  element  techniques  to  model  thermal  oxidation  of 
silicon  in  two  dimensions.  At  temperatures  encountered  in  thermal  oxidation,  sili¬ 
con  dioxide  fiows  viscoelastically.  A  reduced-dimension,  generalized  boundary  element 
method  for  modeling  such  a  problem  has  been  developed.  With  a  Laplace  transform 
technique,  a  viscoelastic  kernel  function  is  derived  from  Kelvin’s  solution,  which  is  the 
fundamental  solution  to  linear  elasticity.  Constant-velocity  loading  is  chosen  to  oper¬ 
ate  with  a  wide  range  of  stress  relaxation  times.  This  scheme  is  capable  of  replacing 
boimdary  element  methods  developed  for  slow  viscous  flow  and  elastic  deformation. 
The  oxidant  diffusion  problem  is  solved  using  a  standard  potential  method  for  Laplace 
problems.  Generated  by  oxide  growth,  stress  affects  both  oxidant  diffusion  and  oxide 
flow.  In  particular,  it  changes  the  diSusivity  of  oxidants  emd  viscosity  of  oxide,  render¬ 
ing  the  diffusion  and  flow  problems  nonhomogeneous.  Domain  solutions  axe  needed  for 
both  processes.  A  unified  initial  stress /built-in  field  formulation  has  been  developed 
to  account  for  the  nonlinear  effects,  using  interior  cells  that  axe  placed  where  stress 
is  significant.  The  interior  solutions  are  realized  with  an  interfacial  source  method, 
whereby  an  area  integral  for  a  cell  is  transformed  to  a  line  intgral  on  the  perimeter 
of  the  cell.  It  has  been  found  that  kernel  functions  based  on  Kelvin’s  solution  are  de¬ 
ficient  in  modeling  incompresible  materials  with  a  “hole”.  A  correction  method  that 
uses  a  source  term  placed  at  the  center  of  the  hole  has  been  implemented  to  overcome 
the  numerical  problem. 
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Chapter  1 


Introduction 


Thermal  oxidation  is  a  principal  process  in  the  fabrication  of  silicon  integrated  circuits 
(ICs).  In  this  process,  silicon  wafers  are  exposed  to  oxidizing  gases  at  elevated  tem¬ 
peratures  to  produce  silicon  dioxide  on  the  surface  of  the  wafers.  Because  of  its  many 
desirable  properties,  such  as  low  surface  states,  electriczd  insulation,  and  masking  abil¬ 
ity  against  certain  dopants,  thermal  oxide  is  used  in  many  applications,  ranging  from 
gate  insulation  for  metal-oxide-semiconductor  field-effect  transistors  (MOSFETs)  to 
dielectric  isolation  between  devices.  Together  with  some  other  processing  characteris¬ 
tics,  thermal  oxidation  currently  makes  silicon  the  most  widely-used  material  for  solid 
state  electronic  devices. 


1.1  Survey  of  Oxidation  Techniques 


Oxidation  on  a  silicon  wafer  can  be  localized  by  patterning  a  layer  of  masking  material 
such  as  silicon  nitride  in  regions  where  no  oxide  is  desired.  Most  MOSFET  ICs  are 
produced  with  the  local  oxidation  of  silicon  (LOCOS)  technique.  A  typical  MOSFET 
structure  is  shown  in  Fig.  1.1  (ij.  To  avoid  excessive  stress  exerted  on  the  silicon 
substrate  during  the  field  oxidation  process,  a  thin  layer  of  pad  oxide  (relief  oxide) 
is  typically  sandwiched  between  the  silicon  nitride  mask  and  the  substrate.  Unfortu- 
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Figure  1.1:  Cross  Section  of  a  MOSFET  Structure.  Only  source  and  gate  of  the  device 
is  shown.  The  field  oxide  region  that  becomes  narrower  as  it  gets  closer  to  the  source 
is  commonly  referred  to  as  the  bird’s  beak. 

nately,  the  pad  oxide  also  promotes  lateral  oxidation  at  the  edge  of  the  nitride  mask. 
As  a  result,  the  termination  region  of  the  field  oxide,  commonly  referred  to  as  the  bird’s 
beak,  is  widened.  Due  to  physical  constraints  atnd  processing  considerations,  the  bird’s 
beak  may  not  be  made  arbitrarily  narrow  without  introducing  serious  side-effects. 

In  order  to  increase  IC  system  complexity  and  functionality,  devices  are  scaled 
down  to  improve  packing  density  and  operation  speed.  Consequently  the  bird’s  beak 
becomes  a  critical  issue  as  it  takes  up  valuable  real  estate.  Its  impact  on  device 
characteristics  also  becomes  equally  important  as,  for  a  small  device,  its  encroachment 
into  the  active  area  has  a  large  effect  on  the  local  dopants  concentration  and  on 
the  fringing  electric  field.  Thus  an  accurate  modeling  of  the  formation  of  the  bird’s 
beak,  which  is  two-dimensional  in  geometry,  will  help  both  process  design  and  device 
modeling.  One-dimensional  process  simulation  progrzuns  csmnot  perform  such  a  task; 
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Figure  1.2:  Sidewall  Masked  Isolation  Structure  (SWAMI).  Basic  processing  steps  and 
evolution  of  the  oxide  shape. 


two-dimensional  simulators  must  be  used. 

One  technique  that  does  not  suffer  from  the  formation  of  the  bird’s  beak  is  the 
sidewall  masked  isolation  (SWAMI)  [2|.  To  create  a  sloped  trench  in  the  silicon  sub¬ 
strate,  it  utilized  an  anisotropic  wet  chemical  etch  that  attacks  in  certain  preferred 
crystallographic  directions.  Then,  a  few  more  process  steps  aire  involved  to  deposit 
and  pattern  silicon  nitride  layers  to  cover  all  surfaces  except  the  bottom  of  the  trench. 
After  the  field-oxidation  step,  the  trench  is  filled  with  oxide.  The  key  steps  are  out¬ 
lined  in  Fig.  1.2.  Good  planarity  is  maintained,  as  it  can  be  seen.  A  bird’  beak  is 
formed  near  the  bottom  of  the  trench,  but  it  does  not  intrude  into  the  active  device 
area,  which  is  located  at  the  top. 

Trench  isolation  is  emerging  as  an  important  isolation  technique  for  high  perfor¬ 
mance  circuitry.  This  process  uses  plasma  etching  to  form  trenches  in  the  silicon 
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Figure  1.3:  Trench  Capacitor  Cell.  After  trench  is  created  by  plasma  etching,  a  thermal 
oxide  layer  is  grown.  Then  a  layer  of  dielectric  (such  as  silicon  nitride)  is  deposited. 
The  hole  is  filled  with  polysilicon  to  serve  as  a  capacitor  plate. 

substrate.  Due  to  directional  selectivity  of  this  etching  technique,  a  high  aspect  ratio 
of  trench  depth  to  opening  width  can  be  achieved,  as  show  in  Fig.  1.3  [3].  For  a  small 
substrate  real  estate,  the  trench  wall  can  provide  a  large  surface  area  suitable  for  use 
as  a  storage  capacitor.  After  thermal  oxide  is  grown  and  some  dielectric  material  is 
deposited,  the  trench  b  filled  with  chemical-vapor  deposition  (CVD)  polysilicon  to 
serve  as  one  plate  for  the  capacitor.  The  substrate  acts  as  the  other  plate.  The  trench 
isolation  is  used  extensively  in  high-density  random-access  memory  (RAM)  such  as  one 
megabit  dynamic  RAM.  This  structure  is  also  suitable  as  a  beu’rier  in  device  isolation 
and  is  being  used  in  bipolar  circuits. 
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Dislocations  on 
(111)  planes 


Figure  1.4:  Schematic  of  a  Structure  with  Dislocations.  The  slips  are  formed  on  (ill) 
planes.  The  lower  defecjis  density  on  the  right-hand  side  is  due  to  stress  relief  provided 
by  the  pad  oxide. 

1.2  Stress  Problems  and  Measurements 

In  the  LOCOS  process,  the  function  of  the  pad  oxide  is  to  reduce  (1)  stress  due  to 
thermal  expansion  nusmatch  between  silicon  nitride  mask  and  the  substrate  and  (2) 
stress  generated  during  the  oxidation  process.  When  high  enough,  stress  damages  the 
underlying  silicon  substrate,  auid  severely  degrades  its  electrical  properties.  Defects 
usually  appear  in  the  form  of  plane  slips  on  {ill}  planes,  as  shown  in  Fig.  1.4.  Stress 
also  affects  the  oxidation  process  itself.  The  phenomenon  of  oxide  thinning  has  been 
observed  in  many  structures.  The  reduction  in  oxide  thickness  occtirs  in  regions  where 
stress  is  significant.  Experimental  characterizations  include  gate  oxide  growth  [4]  and 
oxidation  of  step-shaped  silicon  structures  [5]. 

Possibly  the  most  comprehensive  experiments  on  quantifying  stress  effects  on  ox- 
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idation  were  done  by  Kao  [6].  Cylindrical  structures  of  different  radii  were  oxidized 
at  different  temperatures.  Oxide  thickness  reductions  on  both  convex  and  concave 
surfaces  were  measured  and  plotted  against  the  curvatures,  as  shown  in  Fig.  1.5.  Be¬ 
cause  cylindrical  geometry  is  simple  and  smooth,  the  oxide  growth  behaviors  are  more 
disposed  to  numerical  analysis  than  other  approaches.  Kao  showed  that  a  model  based 
on  incompressible  creeping  flow  of  oxide  can  be  fitted  to  his  experimental  data,  for 
temperatures  as  low  as  800®C[7]. 

Although  many  processing  techniques  have  been  tried  to  reduce  stress  level,  one 
thing  has  yet  to  be  accomplished  is  the  measurement  of  actual  stress  in  local  oxida* 
tions.  Intrinsic  stress,  which  is  due  to  volumetric  expansion  of  oxidized  silicon,  has 
been  deduced  by  measuring  the  degree  of  bending  of  silicon  wafers  m  situ  [8].  The 
viscoelastic  relaxation  of  stress  can  be  detected  with  X-ray  topography  that  has  a 
resolution  as  fine  as  1  p,m  [9].  However,  it  does  not  provide  stress  values.  Raman 
scattering  [10,11]  and  electron  channeling  pattern  [12]  are  used  to  measure  residual 
stress  in  the  silicon  lattice,  after  oxidation.  The  level  of  stress  generated  during  local 
oxidations  is  still  a  matter  of  speculation.  This  is  problematic  because  stress  values 
produced  by  numerical  simulations  are  scalable,  as  explained  in  Chapter  5. 

1.3  Thesis  Objective 


The  goal  of  undertaking  this  thesis  is  to  advance  numericzd  techniques  for  modeling 
thermal  oxidation  in  two  dimensions.  The  boundary  element  method  (BEM),  which  is 
also  known  as  the  boundary-  integral  equation  method  (BIEM),  deals  with  an  integral 
form  of  the  governing  equations.  It  is  particularly  attractive  for  certain  classes  of  linear 
physical  problems  becaxise  it  offers  reduced  dimensionality  in  the  system  solutions. 
Only  boimdaries  are  segmented,  the  enclosed  domain  does  not  need  any  partitioning 
as  in  the  finite  element  method  (FEM).  Hence  the  simulation  setup  process  is  simple. 


Tasks  of  the  thesis  include  the  extension  of  the  incompressible  viscous  flow  model. 
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Figure  1.5:  Experimental  Data  for  Cylindrical  Structures,  (a)  A  cylindrical  silicon 
ring,  (b)  Top  view  of  the  ring.  After  oxidation,  oxide  thickness  is  measured  on  both 
the  concave  and  convex  surfaces  of  the  cylinder,  (c)  The  vertical  axis  is  the  oxide 
thickness  normalized  to  that  of  a  flat  surface.  The  horizontal  axis  is  the  inverse  of  the 
radius  of  the  silicon  boundary. 
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developed  previously  by  the  author  for  his  Master’s  thesis,  to  include  models  such 
as  elastic  deformation  and  viscoelastic  flow.  A  method  based  on  Laplsure  transform 
technique  is  Anally  developed.  It  can  deal  with  such  a  wide  range  of  shear  stress 
relajcation  times  that  it  has  supplanted  the  elastic  and  viscous  flow  formulations  as  the 
default  model  in  this  project. 

Another  objective  is  to  investigate  methods  for  incorporating  nonlinear  effects  in 
the  boimdary  element  method.  Unlike  the  flnite  element  or  flnite  difference  methods, 
the  BEM  does  not  readily  model  nonlinear  problems.  Enhancements  or  approximations 
are  needed  to  deal  with  them.  Several  attempts  had  made  to  approximate  the  nonlinear 
behaviors  so  as  to  avoid  domain  calculations,  but  they  did  not  provided  satisfactory 
results.  Finally,  a  unifled  method  utilizing  interior  cells  is  developed  for  both  the 
diffusion  and  viscoelastic  flow  problems.  The  reoige  of  nonlinearity  it  can  handle  is 
adequate  for  oxidation  modeling. 

1.3.1  Modeling  Considerations 

The  dynamics  that  sustain  thermal  oxidation  zure  (1)  the  supply  of  oxidizing  species 
to  react  with  silicon  to  form  silicon  dioxide,  and  (2)  the  motion  of  the  newly  created 
oxide.  The  main  transport  mechanism  for  oxid2uit  b  diffusion  by  which  the  oxidant 
moves  from  the  free  oxide  sxirface,  through  the  oxide,  to  the  oxide-silicon  interface 
where  it  reacts  with  silicon.  Due  to  volume  expansion,  newly  formed  oxide  partly  fllb 
the  void  left  behind  by  consumed  silicon,  and  peirtly  displaces  exbting  oxide  towards 
the  surface.  The  role  of  masking  materiab  like  silicon  nitride  b  to  block  the  supply  of 
oxidants,  thus  halting  the  oxidation  process  locally. 

A  simple  one-dimensional  oxidation  model  has  been  around  for  about  twenty  years. 
Developed  by  Deal  and  Grove  [13],  thb  b  a  diffusion  model  that  considers  only  the 
transport  of  oxidants,  while  oimtting  that  of  the  oxide.  Its  accuracy  in  predicting 
oxide  thicknesses  has  been  verified  for  many  oxidation  conditions.  In  one-dimensional 

■  tf 
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problems,  mechanics  of  oxide  motion  is  not  considered  because  it  is  uniform  everywhere 
and  therefore  can  be  treated  as  a  translation  in  the  space  dimensions.  This  uniformity 
condition  is  violated  in  the  bird’s  beak  region,  or  in  any  place  that  is  not  planar.  Under 
such  circumstances,  one  must  consider  the  non-uniform  motion  of  oxide.  Viscous  flow 
or  viscoelastic  flow  are  two  models  that  are  commonly  used  for  describing  the  oxide 
motion. 

Research  efforts  in  modeling  thermal  oxidation  of  silicon  can  be  divided  in  two 
camps  according  to  the  numerical  technique  -  the  finite  element  method  and  the 
boimdary  element  method.  For  the  oxidant  diffusion  problem,  all  formulations  zire 
based  on  Deal-Grove  model.  For  the  oxide  motion  problem,  all  three  possible  models, 
namely  elastic  deformation,  viscoelastic  flow,  and  viscous  flow,  have  been  considered. 
The  finite  element  efforts  concentrate  on  viscous  flow  model  whereas  the  boundary 
element  works  emphasize  on  viscoelastic  flow.  Finite  element  approaches  will  be  re¬ 
viewed  first,  followed  by  the  BEM.  A  survey  of  applications  of  the  BEM  method  on 
related  disciplines  will  be  given  at  the  end. 

1.3.2  History  of  Finite  Element  Approaches 

The  finite  element  method  (FEM)  Is  used  widely  to  analyze  mechzLnical  structures, 
fluid  flows,  and  other  continua.  Thermal  oxidation  poses  a  special  problem  for  the 
FEM  in  that  the  oxide  bulk  deforms  and  expands  continuously.  The  mesh  that  divides 
the  bulk  into  elements  must  be  generated  successively  with  an  automatic  algorithm. 
Thus  an  important  task  in  applying  the  FEM  to  thermal  oxidation  is  to  develop  a 
robxist  mesh  generator. 

In  using  the  FEM  to  simulate  local  oxidation,  Poncet  considered  a  linear  elastic 
model  and  a  lairge-strzdn  nonlinear  viscoelastic  flow  model  to  calculate  oxide  motion 
[14].  Unfortunately,  no  stress  values  were  provided  to  judge  the  nximerical  aspect  of 
his  approach.  Poncet  had  also  investigated  the  effect  of  "mechanical  stress  potential” 
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on  diffusion.  Rafferty  incorporates  the  incompressible  viscous  model  for  the  oxidation 
part  of  a  general-purpose  process  simulator,  SUPREM-IV  [15|.  Sutaxdja  treated  oxide 
as  an  incompressible  nonlinear  viscous  fluid  in  the  oxidation/reflow  simulator  CREEP 
[16].  Sutardja  also  proposed  modifications  to  Kao’s  model  for  modeling  stress  effects 
on  oxidation. 

1.3.3  Boundary  Element  Formulations 

The  bo\mdary  element  method  has  gained  the  acceptance  as  a  tool  for  studying  struc¬ 
tural  problems.  Its  greatest  claim  to  fame  is  its  reduced  dimensionality  feature  on 
modeling  linear  problems.  Because  it  only  requires  the  outline  definition  of  a  struc¬ 
ture,  the  analysis  process  is  greatly  simplified,  especially  for  three-dimensional  prob¬ 
lems.  However,  a  major  drawback  of  the  BEM  is  that  it  is  primarily  intended  for  linear 
systems.  Techniques  for  dealing  with  nonlinearity  are  not  fully  developed,  and  their 
efiiciencies  have  yet  to  be  fully  demonstrated. 

The  first  general-ptirpose  two-dimensional  oxidation  simulator  was  developed  by 
Chin  [17|.  Chin  used  the  boundary  element  method  to  model  both  diffusion  and 
flow.  However,  because  the  viscous  flow  problem  was  solved  as  a  coupled  Poisson 
systems,  domain  computations  were  needed.  Later,  the  author  of  this  report  showed 
that  viscous  incompressible  flow  is  a  biharmonic  system  that  is  subject  to  the  same 
reduced-dimension  treatment  as  the  diffusion  problem  [18].  Matsumoto  investigated 
a  viscoelastic  BEM  that,  like  Chin’s  approach,  needs  a  mesh  to  subdivide  the  domain 
[19].  The  bending  of  the  silicon  nitride  mask  was  modeled  as  zm  elastic  process.  But 
despite  the  large  deformation  of  the  thick  nitride  mask,  reported  stress  values  were 
relatively  small.  This  suggests  the  possibility  that  the  stress  history  is  not  updated 
correctly.  Although  its  constraints  are  less  stringent  than  those  of  a  FEM  mesh,  a 
dense  BEM  mesh,  as  \ised  by  Chin  and  Matsiimoto,  severely  degrades  the  speed  per¬ 
formance.  To  avoid  domain  calculations,  Isomae  used  a  Laplace  transform  technique 
to  model  viscoelastic  flow  of  oxide  [20].  However,  numericzJ  approximations  to  the 
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inverse  transform  are  needed. 

1.3.4  Other  BEM  Developments 

Thermal  oxidation  is  a  rather  limited  field  for  the  BEM.  In  other  areas,  many  advances 
have  been  made.  The  following  are  relevant  to  our  resezirch  efforts.  Kaneko  consid¬ 
ered  the  application  of  the  Laplace  transform  technique  to  modify  Kelvin’s  solution 
(21]  for  viscoelastic  deformation,  much  the  same  as  our  approach.  However,  because 
of  the  narrow  interpretation  of  the  constitutive  equations,  the  pzurticular  solution  is 
not  adaptable  to  viscoelastic  flow.  A  new  solution  needs  to  be  derived.  In  a  different 
approach,  Tanaka  uses  the  original  Kelvin’s  solution,  and  makes  the  elastic  modu¬ 
lus  and  Poisson’s  ratio  time-dependent  functions  [22].  For  nonhomogeneous  Laplace 
problems,  Butterfield  suggested  a  “direct  method”  technique  [23].  Om  method  for 
nonhomogeneous  Laplace  problems,  which  is  modeled  after  elastoplastic  techniques, 
turns  out  to  be  similar  to  Butterworth’s,  except  that  ours  is  an  “indirect”  formulation. 


1.4  Report  Outline 


Other  than  the  fact  that  oxide  motion  is  a  “higher”  order  system  compeired  to  oxidant 
diffusion,  these  processes  share  many  common  features.  Therefore  they  will  be  pre¬ 
sented  in  parallel,  from  the  physics,  the  mathematics,  to  the  numerical  approximation. 

Chapter  2  describes  the  basics  of  thermal  oxidation.  The  oxidant  diffusion  and 
oxide  motion  issues  will  be  dealt  with  in  order;  the  physics  will  be  described,  but  the 
mathematics  will  not  be  covered  extensively.  Fundamental  concepts,  especially  those 
of  solids,  will  be  introduced.  Modeling  aspects  of  stress  effects  will  be  discussed. 

The  mathematical  backgrounds  are  extended  in  Chapter  3  to  bring  out  the  charac¬ 
teristics  of  the  systems.  Their  relation  to  the  boundary  element  will  be  studied.  The 
foundation  of  Laplace  and  bihaumonic  problems  will  be  covered.  The  use  of  Laplace 
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tramsform  on  deriving  the  viscoelastic  kernels  will  be  illustrated.  A  unified  approach 
on  modeling  nonlinearities  in  diffusion  and  viscoelastic  fiow  will  be  presented  at  the 
end  of  the  chapter. 

Chapter  4  is  on  the  numerical  realization  of  the  integral  equations.  Basics  such  as 
boimdary  and  domain  segmentation,  methods  for  evaluating  line  integrals,  and  matrix 
solution  for  multiple  domains  will  be  covered.  Adaptive  iterative  techniques  required 
for  solving  nonlinear  behaviors  will  also  be  treated. 

Simulations  of  various  problems  are  shown  in  Chapter  5.  Stress  behavior  is  studied 
as  a  function  of  rel2Lxation  time.  Parametric  fittings  on  Kao’s  experimental  data  will 
be  presented,  along  with  a  discussion  on  the  scalability  of  psirameter  values.  Effects 
of  stress  on  viscosity  of  oxide  and  diffusivity  of  oxidants  will  be  demonstrated. 

Finally,  Chapter  6  summarizes  work  that  has  been  performed  and  proposes  possible 
improvements  and  directions  for  future  investigations. 


Chapter  2 


Thermal  Oxidation  Process 


Thermal  oxidation  is  typically  carried  out  at  800-1100  ®C.  Silicon  wafers  are  exposed 
to  an  oxidhing  ambient  such  as  dry  oxygen,  “wet”  oxygen,  or  steam.  Wet  oxygen  is 
prepared  by  passing  oxygen  through  heated  water  to  saturate  it  with  water  vapor  or 
by  direct  burning  of  hydrogen  gas  in  an  oxygen  ambient.  From  the  following  chemical 
reactions,  silicon  dioxide  is  formed: 

O2  +  Si  —*  SiOj 
2HjO  +  Si  — ►  Si02  +  2H2. 

Because  the  oxidation  rate  of  H2O  is  much  higher  than  O2,  thick  oxides  are  typically 
grown  in  wet  oxygen  and  steam  ambients.  Dry  oxygen  is  used  for  growing  thin  layers 
such  as  gate  insulator,  where  good  oxide  quality  is  of  utmost  importance. 

In  a  one-dimensional  thermal  oxidation  model,  the  processes  of  oxidant  diffusion 
and  oxide  movement  are  often  lumped  together  to  form  a  general  expression  that  re¬ 
lates  oxide  thickness  to  oxidation  time.  Extension  of  such  an  approach  to  two  or  three 
dimensions  is  not  feasible.  The  direction  of  oxidant  diffusion  has  no  bearing  on  the  ox¬ 
ide  movement,  therefore  the  two  processes  miist  be  considered  separately,  even  though 
they  are  tightly  coupled.  Their  mathematical  modeb  will  be  presented  to  serve  as  a 
foundation  for  the  boundary  element  method  described  in  the  next  chapter.  A  review 
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on  the  one-dimensional  oxidation  model  is  given,  followed  by  an  examination  of  the 
two-dimensional  oxidant  diffiision  model.  In  the  second  part  of  this  chapter,  a  brief 
sTirvey  of  the  mechanics  of  solids  and  fluids  will  given.  Elastic  deformation  and  incom¬ 
pressible  viscous  flow  will  be  covered,  followed  by  a  generalization  to  viscoelasticity. 
The  influence  of  stress  on  the  oxidation  process  will  also  be  discussed. 


2.1  Diffusion  Mechanism 

The  most  widely-used  oxidation  model  was  formulated  by  Deal  and  Grove  [13].  This 
model  primarily  concerns  with  the  transport  of  oxidants  from  the  gas  phase  to  the 
oxide  and  within  the  oxide,  and  the  chemical  reaction  that  takes  place  at  the  silicon 
interface.  The  kinetics  of  oxide  motion  is  largely  ignored.  The  Deal-Grove  model 
has  been  found  to  be  generally  valid  for  temperatures  between  700  and  1300  ®C,  in 
atmospheric  or  reduced  pressure  [24] .  In  this  model,  steady-state  molecular  diffusion 
is  assumed  to  be  the  main  transport  mechanism  for  the  oxidizing  species.  The  steady- 
state  approximation  is  justified  by  the  low  solid  solubility  of  oxidants  in  silicon  dioxide. 
As  a  consequence  of  their  low  concentrations,  oxidants  must  move  many  times  faster 
than  the  silicon  interface  in  order  to  maintain  oxidation;  thus  their  relaxation  time 
is  much  smaller  than  that  of  the  moving  boundztries.  To  the  oxidants,  the  silicon 
interface  appears  to  be  stationary;  to  the  silicon  interface,  the  oxidzmt  flux  appears  to 
be  in  steady  state.  A  dimensional  analysis  has  been  performed  by  Chin  to  prove  that 
the  steady-state  approximation  on  the  oxidant  diffusion  process  is  indeed  valid  [ij. 

2.1.1  One-Dimensional  Diffusion 

Fig.  2.1  shows  a  one  dimensional  oxidant  profile.  Oxidizing  species  in  the  gas  phase 
enter  the  oxide  layer  at  the  free  surface;  the  incoming  flux  is  given  by 


Ft  =  A(C-  -  C,) 


(2.1) 
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Figure  2.1:  One-Dimensional  Oxidation  System.  Oxidizing  species  enter  from  the 
gas  phase,  diffuse  through  the  oxide,  and  react  with  silicon  at  the  silicon  interface. 
The  bulk  concentrations  at  the  surface  and  at  the  silicon  interface  axe  C#  and  C,-, 
respectively. 


where  h,  C*,  and  C,  are  the  gas-oxide  mass-transfer  coefficient,  the  equilibrium  bulk 
concentration  (when  there  is  no  diffusion),  and  the  surfawre  concentration,  respectively. 


The  fl\ix  in  the  oxide  is  given  by 


Fi  =  — D— — 


=  D 


dx 
{Co  -  C.) 
do 


(2.2) 


where  D,  C,-,  and  d  are  the  diffusion  coefficient,  the  oxidant  concentration  at  the  silicon 
interface,  and  the  oxide  thickness  respectively. 


JF3  is  the  flux  reaching  the  silicon  interface  where  chemical  reaction  takes  place 
to  form  oxide.  The  quantity  of  incoming  flux  must  equal  the  amount  consumed  by 
the  oxidation  process.  Oxide  growth  is  assumed  to  be  given  by  k,Ci,  where  k,  is  the 
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rate  of  siirface  reaction  and  C,-  is  the  interface  concentration.  We  thus  have  the  third 
relationship: 

Fs  =  k,Ci.  (2.3) 

In  steady  state,  =  Fj  =  jPs  as  fluxes  have  to  be  conserved.  Solving  the  three 
simultaneous  equations,  expressions  for  Ci  ajid  C<,  axe  obtained. 


C.  = 


C.  = 


1  , 


Combining  Eqs.  2,3  and  2.4,  the  following  differential  equation  that  determines  the 
oxide  growth  rate  results 


k,  ,  Kdo 


1  +  -^  +  -^ 

where  N  Lb  the  number  of  oxidant  molecules  required  to  form  a  unit  volume  of  oxide. 
The  solution  of  Eq.  2.6  is  the  well-known  linear-parabolic  oxidation  equation  which 
relates  oxidation  time,  t,  to  oxide  thickness,  dg: 


+  AdQ  =  +  r) 


where 


A  =  2D 


2DC- 

Adi  +  (if 
^  B' 

T  is  the  equivalent  time  required  to  produce  the  existing  oxide  thickness  d,-  on  bare 
silicon.  Usually  h  is  removed  from  the  expression  for  A  because  it  is  much  larger  than 
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Parabolic 

Linear 

BIA  = 

(111)  Silicon 

(100)  Silicon 

Dry  O2  C,  *  7.72  x  !0-  C,  =  3.71  x  lO'-  fim/h 

C2  =  6.23  X  10*  p.m/h  All  other  parameters  same  as 

for  (1 1 1)  silicon 

£,  =  1.23  eV 
£2  =  2.0  eV 


H2O  (640  torr)  C,  =  3.86  x  10^  pim^  C2  =  0.97  x  I0«  jim/h 

C2  =  1.63  X  10*  pim/h  All  other  parameters  same  as 

for  (1 1 1)  silicon 

£,  =  0.78  eV 
£2  *  2.05  eV 


Tabl«»  2.1:  Coefficients  for  Linear  and  Parabolic  Rates. 

There  are  two  limiting  cases  to  the  general  relationship,  Eq.  2.7.  When  the  oxide 
is  thin,  the  growth  rate  is  approximately  constant  since  it  is  limited  by  the  reaction 
rate  at  the  silicon  interface;  do  is  then  given  by 

B, 

do«— (t  +  r).  (2.8) 

As  the  oxide  grows  thicker,  the  growth  rate  becomes  limited  by  the  diffusion  rate  of 
the  oxidant  in  the  oxide.  In  this  case,  the  growth  rate  is  inversely  proportional  to 
the  oxide  thickness;  the  relationship  between  d  and  t  is  then  given  by  the  parabolic 
equation 

dlsi  B(t  +  T).  (2.9) 

Hence  B/ A  and  B  are  known  as  the  linear  and  the  parabolic  oxidation  rate  constants 

respectively.  The  values  oi  Bj A  and  B  in  dry  and  wet  oxygen  are  given  in  Table  2.1 

|25]. 

The  Deal-Grove  model  has  been  found  to  apply  very  well  to  a  wide  range  of  condi¬ 
tions.  However,  it  does  not  account  for  the  rapid  initial  growth  phase  in  dry  oxygen: 
when  oxide  thicknesses  is  less  than  SOOA,  anomalous  high  oxidation  rate  is  observed. 
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Another  behavior  not  predicted  by  the  Deal-Grove  model  is  the  nonlinear  power-law 
dependence  of  the  linear  rate  rate  on  the  ambient  pressiire  in  high-pressure  dry  ox¬ 
idation.  Various  modifications  to  the  theory  have  been  proposed  to  explain  these 
phenomena  [26,27,28]. 

2.1.2  Oxidant  Diffusion  in  Two  Dimensions 

The  extension  of  the  diffusion  model  to  two  and  three  dimensions  is  straight  forward. 
The  flux  vector  given  by  Pick’s  first  law  is  ^ 


F  =  -DVC  (2.10) 

At  the  oxide  free  surface,  Eq.  2.1  is  replaced  by 

T-t  =  -h{C*-C), 


or  equivalently  by 


—  8C  . 

=  h{C  —  C). 


where  n  is  the  unit  normal  vector  at  the  siirface  pointing  away  from  the  bulk.  As 
mentioned  earlier,  the  mass-transport  coefficient  h  is  large;  hence  C  is  ^proximately 
C*  at  the  free  surface.  For  all  practical  pxirposes,  Eq.  (2.8)  can  be  replaced  by 

C  =  C*. 

At  the  silicon  interface,  Eq.  2.3  is  generalized  to 

F-fi  =  fc,C 


or 


^  n  "  —  kgC. 

on 


(2.11) 


^FVom  now  on,  we  do  not  use  subscripts  on  C  to  denote  the  boundary  type.  This  change  is  to  be 
consistent  with  other  expressions. 
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To  complete  the  model,  the  boimdary  condition  at  the  silicon  nitride  interface 
is  considered.  Being  highly  refractory  and  impermeable  to  oxidants,  silicon  nitride  is 
used  extensively  as  a  masking  material  for  local  oxidation.  Since  no  significant  amount 
of  oxidants  diffuse  through  the  nitride  layer,  the  normal  component  of  the  oxidant  fiux 
at  the  nitride  interface  is  assumed  to  be  zero: 


or  simply 


F  •  n  =  0, 


(2.12) 


2.1.3  Ejffects  of  Stress  on  Oxidant  Transport 

Stress  affects  the  diffusion  of  oxidants  by  way  of  several  mechzmisms.  The  ones  that 
Kao  implemented  in  his  oxidation  model  are  on  the  reaction  rate  and  the  equilibrium 
concentration  [29].  According  to  Kao’s  theory,  the  normal  surface  traction  (i.e.  the 
force  acting  perpendicular  to  the  surface)  changes  the  reaction  kinetics  at  the  silicon 
interface.  A  positive  force  impedes  the  flow  of  oxide.  Since  extra  energy  is  required 
to  overcome  the  force,  the  likelihood  of  an  oxidation  event  for  this  energy-activated 
process  is  reduced.  Similarly,  hydrostatic  pressure  can  alter  the  equilibrium  concentra¬ 
tion,  C*,  by  reducing  the  concentration  when  it  is  compressive.  Hydrostatic  pressure 
is  known  to  interfere  with  diffusivity  because  it  removes  vatcancies  needed  for  diffusion 
[30|.  Although  described  but  not  used  in  Kao’s  original  model,  the  diffusivity  model 
was  incorporated  in  Sutairdja’s  simulation  model  [I6j.  Finally,  it  should  be  noted  that 
an  elastic  potential  gradient  can  drive  oxidants  from  from  a  high  pressure  area  to  a 
low  pressure  area  [30,14];  this  process  is  related  to  that  on  C*. 

Followed  Sutardja’s  model,  only  effects  on  fc,  and  D  are  considered  in  the  stress 
study.  Although  other  mechanisms  may  be  included,  they  only  make  the  situation 
worse  by  introducing  more  unknown  variables  than  it  is  currently  possible  to  resolve. 
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The  stress  formulae  are  given  below; 

k,  =  ksoexpi^^^)  (2.13) 

D  =  Doexp(^^)  (2.14) 

where  fc,  and  Dq  are  the  original  unstressed  reaction  rate  and  diffusion  coefficient, 
p„  is  the  normal  component  of  the  surface  traction,  is  the  pressure,  Vjt,  and  Vj 
axe  the  activation  volumes  for  the  reaction  rate  zmd  diffusivity.  The  equations  axe 
given  in  the  activation-energy  form,  to  be  consistent  with  other  temperature-dependent 
expressions.  Note  that  k^  and  Dq  are  functions  of  temperature  as  well. 


2.2  Oxide  Motion 


In  oxidation,  the  silicon  material  undergoes  a  phase  change,  from  single-crystal  (or 
sometimes  polycrystalline)  silicon  into  an  amorphous  glass  which  has  an  open  network 
structure  and  which  does  not  show  any  long-range  crystal  order.  When  a  silicon  atom 
is  oxidized  to  form  a  silicon  dioxide  molecule,  its  volume  increases  from  20  A®  to  45 
A®;  the  expansion  is  approximately  2.25  times.  Part  of  the  oxide  molecule  fills  in  the 
void  of  the  silicon  atom,  the  remaining  squeezes  into  the  existing  oxide  network.  On 
an  atomic  level,  a  high  stress  is  produced.  This  intrinsic  stress  is  responsible  for  the 
bending  of  silicon  substrate  under  planar  oxidation  conditions  [3I|. 

As  oxidation  progresses,  the  silicon  interface  recedes  into  the  silicon  substrate,  while 
the  free  surface  moves  away  from  its  original  position,  eis  shown  in  Fig.  2.2.  From  a 
macroscopic  view,  the  motion  of  oxide  in  one  dimension  can  be  treated  as  a  coordinate 
translation  that  does  not  produce  stress  and  therefore  does  not  require  any  knowledge 
about  the  mechanical  properties  of  oxide. 

For  nonplanax  surfaces  and  nonuniform  oxidation  rate,  the  mechanical  properties 
of  oxide  axe  needed  to  define  a  model  for  the  oxide  motion.  Interactions  with  other 
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Figure  2.2:  Growth  of  Silicon  Dioxide. 

materials  present  in  the  system,  namely  the  silicon  substrate  and  the  silicon  nitride 
mask,  mtist  also  be  contemplated  in  two  and  three  dimensional  models. 

2.2.1  Material  Properties 

The  two  solid  materials  that  oxide  comes  in  contact  are  the  silicon  substrate  and  the 
nitride  mask.  The  properties  of  these  three  materials  are  dissimilar  and  may  require 
different  numerical  treatments. 

Silicon  dioxide  has  been  studied  extensively  for  its  mechanical  properties  [32,33]. 
(But  most  of  these  properties  axe  not  relevent  to  thermal  oxidation.).  Silicon  dioxide 
is  also  known  as  silica;  it  is  usually  found  in  amorphous  state.  At  high  temperatures, 
it  has  a  tendency  to  devitrify  into  crystalline  form.  Like  any  other  materials,  silica 
exhibits  nonlinear  behavior  in  its  elastic  properties.  The  bulk  modulus  changes  with 
applied  pressure.  Under  suitably  high  pressure,  silica  can  be  permanently  densified. 
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Temperature  in  decrees  Centigrade 


1400  1300  1200  !IOO  1000 


Table  2.2;  Equilibrium  viscosities  for  various  silica  glasses.  (A)  type  I  (IR);  (x)  type  II 
(OG  with  0.027wt%  OH);  (•)  type  H  (OG  with  0.04wt%  OH);  (□)  type  HI  (Spectrosil) 
silica  glass. 

This  phenomenon  is  also  observed  in  low-temperature  oxidation  [34]. 

Oxide  is  a  viscoelastic  material  that  flows  readily  at  high  temperatures.  Its  viscosity 
is  a  strong  function  of  the  temperature.  Shown  in  tables  2.2  are  viscosities  for  silica 
glasses  with  different  OH  concentrations  [33,35].  As  it  is  evident  from  the  table, 
viscosity  is  also  very  sensitive  to  the  presence  of  water  vapor.  Water  molecules  combine 
with  bridge  oxygen  ions  in  the  oxide  to  form  stable  nonbridging  hydroxyl  groups  [25], 
thereby  weakening  the  network  structure  and  greatly  reducing  the  viscosity.  Impurities 
such  as  Na  and  P  show  similar  effects.  However,  as  mentioned  in  Chapter  1,  oxide  is 
frequently  approximated  as  a  viscous  incompressible  fluid  for  simplicity. 

Silicon  nitride  is  usually  deposited  by  means  of  chemical  vapor  deposition.  The 
chemical  structure  of  pure  silicon  nitride  is  SijN^.  Like  oxide,  nitride  is  amorphous. 
Its  high-temperature  behavior  is  not  well  characterized,  although  it  is  known  to  crack 
readily.  In  LOCOS,  the  nitride  mask  bends  during  the  process.  When  the  oxide  is 
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removed  by  etching,  the  mask  only  returns  partially  to  the  original  position.  It  is 
suspected  that  silicon  nitride  deforms  elastoplastically  (plastically).  As  a  first  cut, 
nitride  is  modeled  as  a  linear  elastic  element.  If  more  information  is  available,  it  may 
be  treated  as  a  nonlinear  viscoelastic  or  plastic  material. 

The  silicon  substrate  on  which  oxide  is  grown  is  a  single-crystal  material.  For  obvi¬ 
ous  reasons,  we  would  not  like  to  allow  the  substrate  to  “flow” ,  deform,  or  otherwise  be 
damaged  structurally.  Defects  in  the  crystal  structure  is  detrimental  to  the  electrical 
characteristics  of  devices  built  on  it.  Silicon  is  therefore  treated  as  an  elastic  material. 
However,  having  a  long-range  and  well-defined  atomic  arrangement  in  the  structure, 
single-crystal  silicon  exhibits  anisotropy  in  its  mechanical  properties  —  different  crys¬ 
tallographic  planes  have  different  elastic  modulus  and  Poisson’s  ratio.  For  simplicity, 
silicon  may  be  treated  as  if  it  is  an  isotropic  material. 

At  this  point,  it  is  clear  that  a  general-purpose  viscoelastic  method  that  can  han¬ 
dle  a  wide  range  of  stress  relaxation  is  highly  desirable.  As  we  will  see  later  in  the 
next  chapter,  viscous  incompressible  flow  and  elastic  deformation  can  be  considered 
special  cases  of  viscoelasticity  in  which  the  relaxation  time  is  small  or  oo  respectively. 
Nonetheless,  all  the  three  models  will  be  examined  in  the  following  sections  to  show 
their  interrelationship. 


2.2.2  Mechanics  of  Solids  and  Fluids 

For  many  physical  systems,  it  is  often  adequate  to  consider  only  two  (or  even  one) 
dimensional  models  to  imderstand  their  fundamentals  eind  to  perform  analysis  work. 
Their  generalization  to  two  or  three  dimensions  is  straight  forward,  as  we  have  seen 
earlier  in  the  diffusion  process.  For  elasticity,  the  converse  is  true.  One  has  to  start 
with  three-dimensional  concepts,  and  work  the  way  down  to  one  or  two  dimensions.  All 
one  or  two-dimensional  problems  are  based  on  a  particular  simplified  three-dimensional 
behavior. 


CHAPTER  2.  THERMAL  OXIDATION  PROCESS 


33 


Mathematically,  diSiision  is  simpler  than  elasticity.  The  former  deals  with  scalars 
and  vectors,  whereas  the  latter  works  with  vectors  and  tensors.  Vectors  specify  dis¬ 
placement  fields  and  forces;  tensors  defines  the  strain  and  stress  quantities  in  the  ma¬ 
terials.  Because  of  the  increased  “rank”  of  parameters,  indicial  notation  is  employed 
in  this  report  to  represent  components  of  interest  and  to  define  repeated  operations. 
When  an  index  letter  appears  only  once  in  a  term,  it  implicitly  takes  on  a  value  range 
of  1  to  JV  where  JV  is  2  or  3  —  the  dimension  of  the  system.  Thus 


Of  =  iihi 


is  a  compact  representation  for  a  set  of  equations: 


ttj  = 

Ci  =  libj. 


If  an  index  letter  appears  twice  in  a  term,  it  means  sunomation  over  the  range  of  1  to 
N.  Hence 


is  imderstood  to  be 


flt  — 


0*  =  6,(cii  +  C22  +  C35) 
repeated  for  all  possible  values  of  t. 


To  be  compatible  with  the  indicia!  notation,  the  principal  axis  of  the  global  Carte¬ 
sian  coordinate,  z,y,  and  z  are  also  referred  to  as  Zi,Z2,  and  Zj  so  that  they  can  be 
indexed  appropriately.  Subscripts  and  superscripts  are  used  throughout  the  report  for 
different  meanings  and  purposes.  In  order  to  avoid  excessive  confusion,  an  index  letter 
(t,y,  or  k)  appears  only  as  a  subscript  and  only  in  lower  case.  The  indicial  notation, 
although  concise,  is  inconvenient  at  times;  hence,  the  vector  notation  will  be  retained 
wherever  it  is  more  suitable.  Note  that  f  and  /,-  sure  equivalent. 
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Figure  2.3:  Stress  Tensor  Diagram. 

2.2.3  Stress  and  Strain  Tensors 

When  a  force  is  applied  to  a  material,  the  material  deforms  or  rezwrts  accordingly.  Like 
the  diffusion  fliix,  the  force  is  a  vector  that  has  a  magnitude  and  a  direction.  The 
resulting  stress  inside  the  material  is  a  second  order  tensor,  however.  A  stress  tensor 
or  dyadic  has  9  components; 

(7ii  012  ffi3 

<T,J  =  <721  <^22  <^23  (2-15) 

<731  <732  <733 

Their  spatial  orientations  axe  illustrated  in  Fig.  2.3.  an, <722,  and  033  acts  in  the 
direction  perpendicular  to  the  principal  planes  (whose  normals  are  xi,X2,  and  X3); 
they  sure  called  the  normal  stress.  ai2,ai3,a2i,a23,a3i,  and  a32  are  shear  stress  acting 
tangent  to  the  planes.  In  summary,  a,-,-  is  a  force  acting  in  the  j  direction  on  plane  t. 
Due  to  moment  conservation,  the  stress  tensor  is  a  symmetric  matrix,  i.e.  a,-,-  =  a,-,-. 
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The  forces  on  the  3  principal  planes  are  therefore  given  by 

ti  =  <TnXi  +  ciaXj  4-  (TiaXi 
ta  =  o'ai^i  +  "b  ^2z^z 

ts  =  <T3iXi  +  csaXa  +  023X3. 


For  a  plane  of  arbitrary  orientation,  the  force  acting  on  it  is  a  linear  combination  of 
the  3  t’s: 


p  =  Jiiti  +  nata  +  nsts 


which  is,  in  the  indicial  notation. 


Pi 


— - 


(2.16) 


The  n’s  are  the  “direction  cosines”  associated  with  the  plane.  n<  is  defined  to  be  the 
cosine  of  the  angle  between  the  nomaal  of  the  plane  and  the  axis  Xi.  Eq.  2.16  is  one  of 
the  fundamental  expressions  in  mechanics  which  relates  surface  traction  (t.e.  surface 
force]  to  the  stress  tensor  of  the  continuum. 


As  diffusion  is  subject  to  flux  conservation,  a  mechanics  problem  obeys  force  con¬ 
servation,  which  is  given  by  the  equilibrium  condition: 


doii 

where  b  is  the  body  force,  which  may  be  due  to  gravity. 


(2.17) 


Finally,  we  have  the  following  definition  for  the  strain  or  deformation  state: 


= 


1  dm  duj 
2^dxj  dxi* 


(2.18) 


where  e  is  the  strain  tensor  and  u  =  uiXi  +  ujX]  is  the  displacement.  This  statement 
reveals  that  stredn  is  a  unitless  parameter  quantifying  the  distortion  of  a  material.  As 
a  linear  equation,  Eq.  2.18  is  for  small  distortions. 
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It  is  useful  at  this  point  to  separate  the  stress  amd  strain  tensors  into  a  spherical 
part  and  and  a  deviatoric  part.  Given  that  the  spherical  (or  mean)  stress  and  strain 
are  rotationally  invariant,  they  are  often  expressed  as  scalars: 


1 

1 

=  gCii, 


(2.19) 

(2.20) 


where  subscript  M  is  capitalized  to  indicate  that  it  is  not  an  index.  The  deviatoric 
stress  and  strain  tensors  are  obtained  from  the  regular  expressions: 


(2.21) 

(2.22) 

where  the  definition  for  the  Kroncker  delta  is: 


—  ^ij  ~ 


^ij  =  1,  (*  =  y); 

=  Of  (*Vi)* 

The  main  reason  for  separating  the  strain  and  stress  tensor  into  these  two  components 
is  that  the  spherical  part  is  involved  in  volume  change  of  the  material  whereas  the 
deviatoric  part  is  not.  Their  behavior  are  expected  to  be  fimdamentally  different. 


2.2.4  Maximum  Shear  Stress 

A  vector,  having  a  direction  and  a  magnitude,  b  easy  to  picture  regarding  its  transla¬ 
tion  from  one  coordinate  system  to  another.  However,  the  transformation  for  tensor  is 
more  diflScult  to  visualize.  Consider  a  rod  aligned  in  the  ili  direction  and  under  ten¬ 
sion.  It  only  has  a  nonzero  stress  tensor  component,  namely  an.  With  this  particular 
global  reference,  it  is  not  clesir  that  there  is  any  shear  stress.  In  fact,  maximum  shear 
stress  does  occur  on  any  planes  that  are  inclined  at  45“  with  respect  to  the  axis  of  the 
rod. 


CHAPTER  2.  THERMAL  OXIDATION  PROCESS 


37 


Maximiim  shear  stress  is  often  used  as  a  criterion  for  determining  the  plastic  yield  , 
threshold  and  other  nonlinear  behaviors.  Its  value  may  be  obtained  from  the  principal 
stresses.  For  any  stress  state,  it  is  possible  to  find  a  local  coordinate  frame  such  that 
only  Oj],  and  are  non2ero.  If  these  three  components  are  sorted  and  relabeled 
such  that  <7/  >  <7//  >  <7;//,  the  maximum  shear  stress,  <75  (of  another  frame)  is  given 
by 

~  (2.23) 

Another  parameter  used  for  calculating  nonlinear  effects  is  the  maximum  distortion 
energy  given  by 

Cy  =  -  [(<7/  —  ffliy  +  (o-/  —  (TiiiY  +  (<7//  —  <7///)*]  .  (2-24) 

Both  parameters  have  been  used  in  this  investigation.  The  calculation  for  the  principal 
stresses  is  given  in  the  Appendix. 

Eqs.  2.16,  2.17,  and  2.18  are  some  basic  concerns  that  do  not  depend  on  the  de¬ 
tails  of  the  material  types.  What  decides  whether  a  material  is  fiuid  or  solid  is  the 
constitutive  equations  that  relates  the  stress  tensor  to  the  strain  tensor.  The  three 
basic  property  models,  namely  elastic,  viscous,  and  viscoelastic,  are  examined  in  the 
following  sections. 

2.2.5  Linear  Elasticity 

A  purely  elastic  material  has  perfect  memory.  When  a  force  is  applied  to  it,  it  deforms. 
However,  it  will  returns  to  its  original  shape  after  the  force  is  removed.  Expressed  in 
terms  of  the  deviatoric  and  spherical  components,  the  constitutive  equation  for  linear 
elasticity  is: 


3ij  =  2Ge,y 

(2.25) 

Cm  =  ZKcm 

(2.26) 

\ 
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where  G  and  K  are  the  shear  and  bulk  moduli.  These  moduli  may  be  given  In  terms 
of  Young’s  modulus,  E,  and  Poisson’s  ratio,  u: 

^  “  3(1 -2>/) 

Young’s  modulus  is  more  suitable  for  uniaxial  analysis  than  for  2  and  3  dimensional 


G  = 


(2.27) 

(2.28) 


problems.  Poisson’s  ratio  is  a  measure  of  the  compressibility  of  a  materizJ;  it  ranges 
from  0  to  0.5.  When  t/  is  0,  the  material  is  said  to  be  totally  compressible.  At 
this  point,  the  decomposition  of  stress  and  strain  into  spherical  and  deviatoric  parts 
is  not  necessary  because  2G  =  ZK.  At  the  other  extreme,  the  material  is  totally 
incompressible  when  i/  is  0.5.  Because  the  bulk  modulus  is  infinite,  the  relationship 
between  spherical  stress  and  strain  as  specified  in  Eq.  2.28  breaks  down,  must  be 
zero  at  all  times,  but  <tm  take  on  zmy  finite  value.  In  other  words,  deformation 
cannot  involved  volume  change  in  the  material. 


When  the  body  force,  b,  is  zero,  Eqs.  2.18,  2.17  and  2.29  establish  a  constraint  for 
the  displacement  field  u: 

AiV*u  +  (A  +  ;i)V(V.u)  =0.  (2.32) 

Eq.  2.32  is  known  as  the  Cauchy-Navier  equation. 

Listed  in  the  table  below  axe  the  mechanical  properties  of  silicon  dioxide,  silicon, 
and  silicon  nitride: 
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Material 

Silicon  Dioxide 

Silicon  Nitride 

Silicon 

Young’s  Moduliis  (JS'),  dynes/cm“^ 
Poisson’s  ratio  {u) 

Bulk  Modulus  {K),  dynes/cm“’ 
Shear  Modulus  (G),  dynes/cm"* 

8  X  10“  ~ 

0.194 

4.35  X  10^^ 

3.35  X  10“ 

3.29  X  10“ 
0.266 

2.34  X  10“ 

1.3  X  10“ 

3.69  X  10“ 
0.42 

7.68  X  10“ 
1.3  X  10“ 

Table  2.3:  Elastic  properties  of  oxide,  nitride,  and  silicon. 


2.2.6  Two-Dimensional  Plane  Strain 

In  making  a  two-dimensional  approximation  for  a  three-dimensional  problem,  two 
approaches  may  be  taken.  In  plane  stress,  the  domain  is  assumed  to  be  thin  and  free 
to  expand  in  the  third  dimension,  say  X3,  and  loading  on  the  boundary  does  not  vaxy 
in  the  third  dimension.  This  results  in  the  following  constraints  on  the  stress  and 
strain  tensors: 

=  €is  =  C32  =  631  =  0  (2.33) 

Cii  =  ^23  =  ^32  —  ^31  —  ^33  —  0  (2.34) 

€33  =  +  0’22).  (2.35) 

In  plane  strain,  the  domain  is  assumed  to  be  long  and  uniform  in  the  Z3direction. 
Because  it  cannot  expand  In  £3  direction,  the  following  constraints  are  imposed  on  the 
stress  and  strain  tensors: 

€i3  =  C23  =  €33  =  C32  =  C31  =  0  (2.36) 

<713  =  (723  —  <7z2  —  ^31  =  0  (2.37) 

^33  =  ^{<^11  +  <^22)*  (2.38) 

The  respective  configurations  for  plane  stress  smd  plain  strain  are  illustrated  in  Fig.  2.4. 

For  our  simulation  efiforts,  we  assiime  all  structures  are  long  in  the  third  dimension 
and  therefore  treat  them  as  2-dimensionaI  plane  strain  problems.  The  Cauchy-Navier 
equation  Eq.  2.32  and  other  essential  relationships  remain  unchanged  for  plane  strain. 
Unless  otherwise  noted,  we  will  refer  to  all  parameters  and  functions  as  if  they  are  for 
two-dimensional  systems  from  here  on. 
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(a)  (^) 


Figure  2.4:  Two-Dimensional  Elasticity,  (a)  Plane  stress;  (b)  plane  strain. 

2.2.7  Incompressible  Viscous  Flow 

Viscous  flow  is  one  model  totally  different  from  eleisticity.  It  is  totally  memoryless  - 
any  deformation  is  permanent.  Stress  is  only  sustained  only  the  material  deforms  con¬ 
tinuously,  i.e.  when  it  flows  as  a  fluid.  The  dyneunical  equation  for  slow  incompressible 
viscous  flow  is: 

t;  V*v  =  -  VP  (2.39) 

where  rj  is  the  viscosity,  P  the  hydrostatic  pressure,  and  v  =  dnfdt  is  the  velocity. 
This  parameter  is  the  same  as  the  spherical  stress  in  el2isticity  but  with  opposite  sign: 
P  =  —Cm-  a  positive  cm  implies  tensile  pressure  whereas  a  positive  P  is  compressive. 
The  flow  is  subject  to  the  incompressibility  condition: 

V-v  =  0  (2.40) 

which  is  the  equivalent  to  the  =  0  statement  for  incompressible  elasticity. 
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2.2.8  Viscoelastic  Deformation 


Showing  intermediate  behavior  between  elastic  deformation  and  viscous  flow  is  a  wide 
range  of  viscoelastic  materials.  They  are  further  classifled  as  Maxwell  fluid  or  Voigt 
solid.  Like  an  elastic  material,  the  Voigt  solid  will  rettirn  to  its  original  shape  after 
the  applied  force  is  removed.  It  differs  in  one  aspect  -  the  strain  is  not  only  related 
to  stress  but  also  to  stress  rate.  The  Maxwell  element  behaves  like  a  viscous  fluid  but 
it  retains  some  memory.  If  the  applied  force  is  removed  fast  enough,  it  will  return 
partially  to  its  original  position.  On  the  other  hand,  if  it  is  kept  in  a  distorted  state 
for  an  extended  period,  it  will  lose  all  its  stress  and  remain  in  the  new  shape. 


Formally,  viscoelastic  stress-strain  relationship  is  given  by  a  time-differential  equa¬ 
tion  of  the  form  "" 

Po<^  +  Picr  +  p2Sr  qoe  +  qii  +  q2e  (2.41) 

or  by  the  more  compact  statement: 

Pff^Qe  (2.42) 


where 


a* 


P  = 

»=0 


i=0 


Table  2.2  show  some  fluid  and  solid  models  [36].  In  general,  the  material  is  a  solid 
if  ^0  is  nonzero.  We  note  that  the  elastic  constitutive  equations  Eq.  2.26  is  a  special 
case  of  Eq.  2.42.  The  time-differential  stress-strain  relationship  can  be  converted  to 
different  forms  to  suit  the  applications  -  hereditary  integral,  complex  variable,  Laplace 
transform,  and  so  forth. 


Although  one  may  incline  to  use  a  high-order  differential  equation  to  describe  a 
viscoelastic  behavior  as  completely  as  possible,  such  an  overkill  approach  may  not 
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serve  to  illuminate  the  problem.  In  fact,  there  may  even  be  difficulties  In  gathering 
sufficient  data  to  estimate  the  basic  p’s  and  g’s.  Factors  such  as  nonllne2Lrlty  deserve 
more  attention  In  advanced  approximations. 

For  our  simulation,  we  only  consider  the  simplest  viscoelastic  flow  representation. 
Here,  we  assume  that  the  spherical  component  Is  purely  elastic: 

Ojvf  =  3K€m. 

And  the  devlatorlc  part  Is  a  Maxwell  fluid  governed  by  a  flrst  order  differential  equa¬ 
tion: 

1  dff  a  _d€ 

G  17  dt' 

Its  corresponding  mechanical  model  is  a  linear  spring  In  series  with  a  viscous  dashpot, 
as  depicted  In  Table  2.2.  Note  that  the  spherical  component  must  display  a  solid 
behavior  -  a  fluid  loses  stress  and  allows  indeflnlte  compression,  producing  a  mass 
conservation  problem. 

2.2.9  Nonlinear  Effects  of  Stress 

All  the  models  described  earlier  are  idezdlzed  linear  approximation.  They  may  be  aug¬ 
mented  to  Included  certain  nonlinear  stress  effects.  Consider  a  elastic  problem  whose 
stress-strain  curve  Is  shown  In  Fig.  2.5.  For  small  stress  (or  strain),  the  relationship 
is  linear  from  O  up  to  P;  it  starts  to  depart  from  the  straight  line  and  become  non¬ 
linear  for  larger  stress.  If  the  material  is  able  to  return  to  O  after  the  applied  stress 
is  removed,  It  Is  nonlinear  elasticity  that  we  are  working  with.  If  the  material  has 
experienced  irreversible  slip  on  the  atomic  level,  it  may  only  return  to  the  position  R. 
In  this  case,  we  axe  dealing  with  elastoplatlclty  or  pleistlclty. 

For  our  modeling  efforts,  we  Ignore  the  nonlinear  elastic  moduli  or  plasticity  —  we 
are  only  concerned  with  the  stress  effects  on  viscosity  in  viscous  and  viscoelastic  flow. 
Note  that  the  derivation  of  the  nonlinear  viscoelasticity  BEM  is  based  on  plasticity. 
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Figure  2.5:  Nonlinear  Elasticity  and  Plasticity.  The  material  is  linear  from  O  to  P;  it 
exhibits  a  nonlinear  behavior  from  P  to  Q.  If  it  able  to  return  to  O  when  the  applied 
stress  is  removed,  the  material  is  elastic.  If  it  C2m  only  return  to  R,  it  is  said  to  suffer 
from  permanent  plastic  deformation. 
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In  Kao’s  viscous  flow  model,  the  following  formulation  is  used  to  modify  the  vis¬ 
cosity  [29]: 

J7  =  »7o  exp  (-^^)  •  (2-43) 

However,  as  pointed  out  by  Rafferty  [15],  this  model  can  induce  positive  feedback  for 
certain  geometries.  What  happens  is  that  a  large  negative  (compressive)  stress  catxses 
the  viscosity  to  increase,  this  in  turns  drives  the  stress  more  negative,  until  the  solution 
blows  up^.  Using  the  same  viscous  flow  model,  Sutau’dja  eliminated  the  instability  by 
changing  Eq.  2.43  to  [16] 

V  = 

Becaxise  the  maximum  shear  stress  cs  is  defined  to  be  a  positive  quantity,  the  viscosity 
will  always  be  equal  or  less  than  its  nominial  value.  Thus  the  system  is  stable.  However, 
in  order  to  get  a  good  fit  with  Kao’s  experimental  data,  Sutardja  was  required  to  use 
a  truncated  function  for  the  diffusivity: 

£l  =  DoXimii(«p(2^),Dmax).  (2.45) 

The  value  for  i?max  sometimes  as  low  as  Dq.  To  constraint  D  so  that  it  does  not 
exceed  its  nominal  value  is  not  physical.  If  compressive  stress  cam  reduce  diffusivity, 
tensile  stress  should  enhance  it.  In  Deal-Grove  oxidation  model,  the  retardation  effect 
of  intrinsic  stress  on  diffusivity  is  probably  incorporated  in  the  vaxioTis  parameters. 
When  the  stress  level  is  reduced,  the  actual  diffusivity  cam  be  higher  tham  its  nominal 
value.  Of  course,  there  is  a  linut  as  to  how  much  diffusivity  cam  be  boosted  by  tensile 
stress  before  other  effects,  such  as  fracturing,  set  in. 


riQ 


<rsVo/2kT 

sinh(o-5V’o/2kr) 


(2.44) 


We  propose  the  following  model  for  viscosity: 

f  crsVrjz/kT 

=  j 


smh{<rsV„,/kT) 

which  is  a  combination  of  Eqs.  2.43  and  2.44.  We  retain  the  pressure  dependency  for 

the  following  two  reasons.  First,  in  Kao’s  experiments,  the  main  difference  between  a 

a  viscoelastic  flow  model  is  used  instead,  the  solution  will  always  be  stable  because  the  worst 
that  can  happen  is  that  the  material  becomes  elastic  when  the  viscosity  goes  to  infinity.  Bat  now  it 
suffers  from  a  different  misbehavior:  a  sudden  solidification  from  a  viscoelastic  fluid  to  an  elastic  solid. 
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concave  surface  and  a  convex  one  is  the  sign  of  the  pressure  term,  as  far  as  stress  effects 
are  concerned.  The  pressure  is  tensile  in  the  convex  case,  yet  no  significant  overshoot 
in  oxide  thickness  is  observed  for  any  curvatures.  This  suggests  that  the  viscosity  of 
oxide  is  less  than  its  normal  value  for  the  convex  structures.  Second,  researchers  have 
proposed  that  shear  action  help  to  densify  silica  under  pressure  by  interlocking  the 
Si02  network  [32].  The  same  interlocking  mechanism,  which  requires  the  assistance 
of  compressive  pressxire,  may  make  it  more  difficult  for  Si-0  bonds  to  break,  hence 
increasing  the  viscosity.  Because  of  geometric  considerations,  shear  stress  is  always 
present  in  oxide.  By  choosing  a  suitable  ratio,  we  can  assure  that  Eq.  2.46  is 

always  stable. 


Chapter  3 


Boundary  Element  Formulations 


The  equations  for  diffusion  and  mechanics  are  presented  in  Chapter  2.  In  this  chapter 
we  shall  examine  and  imderstand  their  mathematical  characteristics  and  significance. 
These  properties  serve  as  the  foundation  for  the  boimdary  element  method  (BEM).  The 
development  of  the  boundary  element  technique  will  be  laid  out  in  detail.  Variations 
of  implementation  techniques  will  be  described. 

It  will  be  shown  that  both  oxidant  diffusion  and  oxide  fiow  are  simple  boundary 
value  problems  (BVP’s).  The  definition  of  a  BVP  is  that  the  solution  within  the 
enclosed  domain  is  uniquely  determined  by  the  boundary  conditions.  Mzmy  physical 
systems  are  BVP’s  in  nature.  The  DC  operation  of  a  bipolar  or  field-effect  transistor 
is  one  such  example.  The  terminals  or  boundary  contact  points  are  the  emitter,  base, 
and  collector  for  the  bipolar  device,  or  source,  gate,  and  drain  for  the  field  effect 
transistor.  For  a  given  set  of  terminal  voltages,  there  will  be  some  well-defined  current 
flowing  through  these  terminals.  However,  the  modeling  of  semiconductor  devices  is 
complicated  and  involved  because  the  systems  are  highly  nonlinear.  Linear  BVP’s 
are  easier  to  solved;  in  p2u^icular,  the  so-called  Laplace  and  biharmonic  problems  are 
especially  amenable  to  the  botmdary  element  method. 

As  its  other  more  generic  name  -  the  boundary-integral  equation  method  -  im- 


I 


47 


CHAPTER  3.  BOUNDARY  ELEMENT  FORMULATIONS 


48 


plies,  the  BEM  involves  the  solution  of  some  integral  equations  at  the  boundary.  No 
domain  calculations  are  necessary,  exception  for  nonlinear  problems.  This  reduced  di¬ 
mension  feature,  i.c.  solving  a  three-dimensional  system  as  a  surface  problem  (which 
is  two-dimensional),  and  a  two-dimensional  system  as  a  line  problem  (which  is  one¬ 
dimensional)  is  very  attractive  for  many  applications. 


3.1  Laplace  Problems 

The  solutions  for  many  physical  problems  of  interest,  such  as  electrostatic  potential  in 
insulators,  steady-state  heat  flow  and  steady-state  diffusion,  satisfy  Laplace’s  equation: 

=  0  (3.1) 

where  $  denotes  the  primary  parameter  of  the  system.  For  the  three  examples,  $  is 
the  potential,  the  temperature  distribution,  or  the  concentration  of  diffusing  species, 
respectively.  The  secondary  parameter  of  these  systems  is  the  “flux” ,  a  vector  element 
defined  as 

F  =  (3.2) 

where  F  is  interpreted  to  be  the  electric  field,  thermal  energy  flux,  or  particle  flux, 
and  D  is  the  dielectric  constzmt,  thermal  conductivity,  or  diffusivity. 

The  Laplacian  condition  arises  from  a  conservation  law  on  the  fi-ux 

V  .  F  =  0  (3.3) 

which  states  that  the  divergence  of  F  is  zero,  t.e.  electric  field,  heat  energy,  or  diffusing 
piU'ticles  must  conserve. 

Substituting  Eq.  3.2  into  Eq.  3,3,  the  following  expression  is  obtained, 

VD- =  0,  (3,4) 


This  equation  becomes  Laplace’s  equation  (Eq,  3.1)  if  is  uniform  so  that  VjD  =  0. 
In  other  words,  $  is  Laplacian  if  the  system  is  homogeneous. 
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The  solution  to  a  Laplace  problem  is  uniquely  determined  by  its  boundary  condi¬ 
tions,  which  may  be  given  in  the  following  forms: 

•  $  itself  — a  Dirichlet  boundary  condition. 

•  the  normal  derivative  of  $  —  a  Neumann  boundary  condition. 

3$ 

•  0:$  +  1  where  ot,Pf  and  7  axe  constants.  This  generalized  expression  is 

called  a  Robin  condition. 


A  fimction  f{x,  y)  that  satisfies  Laplace’s  equation  for  two  dimensions  is  also  known 
as  a  harmonic  function,  a  term  ostensibly  comes  from  complex  analysis^.  The  mean 
value  theorem  states  that  the  value  of  a  harmonic  function  at  point  p  is  equal  to  the 
average  of  its  values  over  the  area  of  any  circle  centered  on  p  [37].  From  this  theory,  we 
obtain  the  statement  that  the  maximum  (or  minimum)  of  /(x,  y)  in  the  domain  fl  must 
occur  on  the  enclosing  boundary  F.  This  consequently  guarantees  the  uniqueness  of 
the  solution:  if  two  functions  /i(x,y)  and  /2(x,  y)  are  harmonic  inside  fl  and  are  found 
to  match  along  F,  they  must  be  identical  within  fi:  if  the  difference  /i{x,y)  —  /2(x,  y) 
(which  is  also  harmonic)  is  zero  on  F,  it  has  to  be  zero  everywhere  in  fl,  as  dictated 
by  the  maximum  principle. 

We  have  thus  seen  the  foundation  of  the  BEM  Laplace  problem:  if  w  ■  can  somehow 
find  a  certain  function  that  is  harmonic  inside  fl  and  that  matches  the  prescribed 
condition  on  F,  we  know  that  we  have  obtained  the  solution  for  fl.  Indeed,  we  never 
need  to  figure  out  what  is  happening  inside  fl  during  the  computation  process,  although 
we  can  if  we  want  to. 


The  determination  of  a  solution  that  fits  the  prescribed  boundary  condition  can  be 
accomplished  in  several  ways.  First,  a  polynomial  series  in  x  and  y  may  be  used.  As 


^Th«  periodic  logarithmic  kernel  mentioned  later  in  this  chapter  is  in  fact  obtained  from  a  complex 
fiinction,  as  it  will  be  shown  later. 
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an  illustration,  five  terms  are  shown  below: 

F{x,  y)  =  A  +  Bx  +  Cy  +  Dxy  +  E{x^  -  y^)  (3.5) 

This  particular  harmonic  function  has  five  unknown  coefficients,  therefore  it  can  be 
iised  to  match  “five”  boimdary  conditions.  If  there  are  fewer  boundary  conditions, 
some  terms  in  Eq.  3.5  have  to  be  removed.  Conversely,  if  there  are  more,  higher 
polynomial  terms  such  as  —  xy*  must  be  added.  If  the  problem  is  periodic,  say  in 
the  X  direction,  the  following  trigonometric  series  may  be  used: 

N 

=  5Z[®nCOs(nx)  +  6„sin(nx)]{c„exp(ny)  +  d„exp(-ny)j.  (3.6) 

n=l 

Eqs.  3.5  and  3.6  are  distinctly  different  -  they  do  not  seem  to  be  able  to  yield  the  same 
result.  This  is  because  the  number  of  “boundary  conditions”  is  assumed  to  be  finite. 
In  actuality,  the  boimdary  conditions  vary  continuously  along  F,  therefore  these  series 
should  be  infinite  and  should  produce  identical  results.  In  a  numerical  solution,  only 
a  fihite  number  of  samples  of  the  boundary  conditions  are  used,  in  conjunction  with  a 
finite  series  approximation. 

The  two  techniques  mentioned  above  have  a  common  drawback,  namely,  the  order 
of  the  series  depends  on  the  number  of  sampling  points.  One  would  like  to  use  an 
approximating  function  that  does  not  have  this  dependency  so  that  the  solution  can  be 
generated  in  a  straightforward  manner.  Green’s  function  fulfills  such  a  requirement.  It 
is  a  characteristics  solution  to  the  governing  equation  of  the  system,  Laplace’s  equation 
in  this  case.  All  BEM  formulations  rely  on  Green’s  function  for  solutions. 

The  BEM  may  be  classified  into  two  different  categories  according  to  its  formula¬ 
tions:  the  direct  method  and  the  indirect  method.  The  indirect  method  is  also  known 
as  the  potential,  source,  or  classic  method.  The  direct  method  also  known  as  Green’s 
third  identity  for  Laplace  problems  and  Somigliana’s  formula  for  elasticity  [38]. 

The  direct  method  is  so  called  because  the  unknown  boimdary  parameters  are 
expressed  directly  in  terms  of  the  known  parameters,  namely  the  boundary  conditions. 
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For  the  indirect  method,  an  intermediary  solution  must  be  explicitly  obtained  before 
the  unknown  boundary  parameters  can  be  computed.  On  computational  efficiency, 
both  methods  are  approximately  the  same  -  they  both  require  a  matrix  multiplication 
and  a  matrix  inversion  (solution),  only  the  order  of  these  operations  is  different. 

The  direct  method  is  preferred  by  many  users  because  its  error-weighing  method 
is  similar  to  that  of  the  finite  element  method  and  also  because  of  its  alleged  better 
accuracy.  With  the  absence  of  intermediary  sources,  it  is  readily  interfaciable  with 
the  finite  element  method  or  the  finite  difference  method.  The  indirect  method,  on 
the  other  hand,  is  conceptually  simpler.  Also,  parameters  of  various  types  are  readily 
extractable  or  implementable,  as  it  will  be  shown. 

From  a  physical  standpoint,  the  indirect  method  is  actually  more  “direct”  th«m  the 
direct  method  because  it  models  what  goes  on  physically  at  the  boundary,  such  as  the 
application  of  electrical  charges  or  loading  forces.  Since  it  is  the  method  that  we  use, 
we  shall  begin  with  a  description  of  the  indirect  method. 

3.1.1  Indirect  Method 

Consider  a  two-dimensional  potential  problem.  Suppose  that  the  potential  distribution 
in  the  enclosed  region,  fl,  shown  in  Fig.3.1  is  to  be  modeled  by  distributed  charges  on 
the  boundary,  F,  [39|.  The  potential  generated  by  these  sources  will  be 

^(p)  =  /  p(q)^’(p  -  q)dr  (3.7) 

Jr 

where  $  is  the  potential  at  p,  p  is  the  charge  distribution  zind  ^*(p  —  q)  =  log  [p  —  q|  is 
the  logarithmic  potential  at  p  due  to  a  point  charge  at  q.  <;^*(p-q)  is  Green’s  function 
for  Laplace  problems  because  it  satisfies  Laplace’s  equation  everywhere  except  at  the 
sourc  point: 

VV(q)  =  2T<(q)  (3.8) 

where  5(q)  is  a  Dirac  delta  function.  Eq.  3.7  is  commonly  referred  to  as  Fredholm’s 
integral  equation  of  the  first  kind.  The  function  is  called  a  kernel  or  a  kernel  function 
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^  Positive 

©  Negative 


Figure  3.1:  Source  Distribution  Technique  for  Potential  Problems 
and  p  is  an  unknown  function. 

Let’s  assume  that  a  solution  to  the  source  distribution  p  has  been  determined.  The 
potential  it  generates  will  be  correct  on  P  sind  inside  fl.  Moreover  the  electric  field  can 
be  readily  obtained,  as  shown  below  in  both  vectorial  and  indicial  forms: 

(3.9) 


If,  instead  of  potential,  electric  field  is  specified  at  the  boundary,  we  have  a  Neu¬ 
mann  problem.  As  the  boundary  condition  is  given  by  5$/5n,  the  equation  for  deter- 


E  =  -V$ 

=  -^p(q)v^*{p-q)dr 


where  we  have  omitted  the  dielectric  constant  for  simplicity. 
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mining  the  source  distribution  is 

—  =  -E-n  (3.10) 

=  /^o(p)^<^-(p-q)dr  (3.11) 

Note  that  a  Neumann  problem  is  ill-specified  in  the  sense  that  if  $(p)  is  a  solution  to 
the  potential  distribution,  $(p)-l-c  where  e  is  an  arbitraury  constant  is  also  a  permissible 
solution.  (As  a  homogeneous  solution,  c  does  not  contribute  to  the  boundary  condition 
d^fdn.)  To  guarantee  a  rmique  solution,  an  additional  normalization  equation  on  p 
is  required: 

^/>(q)dr  =  0.  (3.12) 

It  has  the  effect  of  preassigning  a  value  for  c. 

What  has  been  described  is  basically  the  boundary  element  method  of  simple- 
source  distribution  [38,40].  It  is  also  called  the  single-layer  formulation.  For  Dirichlet 
problems,  the  double-layer  method  of  the  following  form  is  sometimes  employed  [38] : 

^(p)  =  ^p(q)^^*io<7(p-q)dr  (3.13) 

=  ^p(q)^^*(p-q)  +  ’rp(p)  (onr)dr  (3.14) 

where  5/3n,  is  the  normal  derivative  at  q  with  respect  to  q.  The  tt  term  in  Eq.  3.14 
results  from  the  integration  over  the  singularity  at  q  =  p.  The  name  “double  layer” 
comes  from  the  observation  that  d<f>* (dn^  is  potential  due  to  a  dipole  charge  whose 
axis  is  normal  to  the  sxirface  -  hence  it  appears  that  the  domain  is  enclosed  by  two 
layers  of  boundary  charges. 

Eq.  3.14  is  known  as  a  Fredholm  integral  equation  of  the  second  kind.  It  produces  a 
diagonally  dominant  matrix  ideal  for  numerical  calculations.  It  can  also  be  solved  very 
efficiently  using  the  Liouville-Neumann  series  [41].  However  the  kernel  functions  be¬ 
come  highly  singular  if  we  attempt  to  compute  V  $  numerically.  Thus  our  application 
uses  the  single  source  formulation. 


CHAPTER  3.  BO  UNDARY  ELEMENT  FORMULATIONS 


54 


Poisson’s  Equation 

The  Laplace  condition  implies  the  domain  is  source-free;  i.e.  there  are  no  electrical 
charges,  heat  sources,  or  diffusing  species  sources  for  the  respective  systems.  If  the 
domain  has  sources,  Poisson’s  equation  applies: 

V*${p)=p{p)  (3.15) 

where  p  is  the  source  density.  Consequently,  Eqs.  3.7  and  3.14  will  have  an  additional 
expression  to  account  for  the  domain  contributions: 

i|^p(q)#-(p-q)dn.  (3.16) 

The  Poisson  condition  often  arises  when  the  domain  is  nonhomogeneou5 .  If  that  is 
the  case,  the  domain  source  p  is  dependent  on  the  boundary  condition. 


3.1.2  Direct  Method 


To  derive  the  direct  method  for  potential  problems,  we  start  with  the  divergence 
theorem: 


/  V.Fdn  =  J^F-Adr. 


(3.17) 


If  we  subtract  the  solution  for  F  =  from  the  solution  for  F  =  where  <j>  and 
V'  are  scalars,  we  obtain  Green’s  second  identity: 


(3.18) 


If  we  replace  ^  and  ip  with  (p*  and  $  defined  earlier  for  the  indirect  method,  we  obtain 
the  direct  formulation: 


—  2?r$(p) 


dr-/^-(q-p)V=4(q)dn, 

Jn 

(3.19) 
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which  is  better  known  as  Green’s  third  identity  or  Green’s  formula.  If  the  domain  is 
source-free,  the  domain  integral  on  the  right-hand  side  of  Eq.  3.19  can  be  removed. 

To  determine  the  unknown  boundary  parameters,  p  is  set  to  T  and  $  is  expressed 
in  terms  of  di/dn  or  vice  versa  *.  One  advantage  of  Green’s  formula  is  that  Eq.  3.19 
does  not  contain  a  source  function.  On  the  other  hand,  it  is  tedious  to  determine 
inside  0  -  it  requires  the  second  derivatives  of  to  do  so.  These  highly  singular 
second  derivatives  are  not  needed  by  the  indirect  method. 

3.1.3  Periodic  Kernels 

So  far,  the  form  of  domain  that  we  are  concerned  with  is  a  simple,  closed  region.  There 
are  many  situations  in  which  one  would  like  to  work  with  a  periodic  structure  that 
extends  from  x  =  — oo  to  x  ==  oo.  In  this  configuration,  the  upper  boundary  Fu  and 
the  lower  boundary  F;  are  separated  and  never  touch.  The  solution  for  $  is  given  by 

^(p)  =  /  Pu(p)<^’(q-p)dF+ /  Pj(P)<^*(q“P)dF  (3.20) 

Vr,,-oo  JTi,-co 

Ass\iming  the  dimensions  are  normalized,  the  following  conditions  will  always  be  sat- 
bfied: 

$(p)  =  $(p  +  2n7r5t) 

P(P)  =  p(p  +  2n7rx) 

for  any  integer  n.  Given  its  periodic  nature,  the  indefinite  integral  can  be  reduced  to 
an  integral  over  a  period: 

f  p(q)<?i*(p-q)dr=  f  p{q)^*{p-q)dF 

OO 

by  introducing  a  periodic  kernel: 

^*(P)  =  '^*(p  +  2n7rx) 

n=— OO 

=  \  log  [2  (cosh(y)  -  cos(x))| , 

description  is  oversimplified.  The  left-hand  side  of  Eq.  3.19  actually  contains  a  scale  factor 
that  depends  on  the  shape  of  the  bonndary.  See  [38,40,42]  for  details. 
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which  is  obtained  from  the  real  part  of  this  complex-variable  expression  [41]; 

lR|log  Jsini(x-f-ty)j|. 

The  singularities  of  are  located  at  2n7r3l  for  any  integer  n.  By  applying  Taylor  series 
expansion  on  cosh(j/)  and  cos(2), 

cosh(y)  sa  H-  y*/2; 
cos(x)  «  1  —  x*/2, 

it  can  be  shown  that  reduces  to  the  non-periodic  version  near  the  origin: 

^(p  «  0)  «  ^’(p). 

For  |y|  :>  0,  we  have  this  approximation: 

^  «  ^|y| 

which  is  Green’s  function  for  one-dimensional  Laplace  problems.  The  far-field  analysis 
is  a  very  useful  tool  because  it  deals  with  a  simplified  one-dimensional  behavior  that 
reveals  the  source  activities  at  y  =  0. 


3.1.4  Source  Formulation  for  Oxidant  Diffusion 


Having  covered  the  principles  of  the  BEM  for  Laplace  problems,  we  now  describe  its 
application  on  modeling  oxidant  diffusion. 


First,  it  will  be  shown  that  diffusion  is  a  Laplace  problem.  Consider  Fick’s  second 
law  of  diffusion: 


dC 

dt 


=  -V-F 


In  steady  state,  dC /dt  =  0.  Thus  the  flux  is  divergence-free: 


V-F  =  -V{DVC) 

=  0. 


(3.21) 

(3.22) 
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If  D  is  uniform  throughout  the  oxide,  as  it  will  be  initially  assumed,  Eq.  (2.12)  is 
reduced  to  Laplace’s  equation: 

V^C  =  0.  (3,23) 


With  the  indirect  BEM  method,  the  oxidant  distribution  is  given  by  this  statement, 

^(P)  =  (3.24) 

where  p  is  interpreted  as  the  ’source’  or  ’sink’  for  oxidants. 

In  the  literature,  parameters  for  oxidation  are  commonly  given  in  terms  of  B  and 
B/A,  not  in  C*,  k,,  and  D.  By  multiplying  all  C’s  with  kjN  and  redefining  p  in 
Eq.  3.24  to  be  pk^fN^  the  oxidant  diffusion  problem  may  be  strictly  expressed  in  terms 
of  B/A  and  B.  Eq.  3,24  now  becomes  a  oxide  ‘growth  rate’  distribution  equation 

where  G(p)  is  the  oxide  growth  rate  if  silicon  is  present  at  p. 


The  boundary  conditions  for  the  free  oxide  surface  and  the  oxide-silicon  interface 
are  replaced  by 

G  =  ^  (3.26) 

A 


AdG 

2‘^  =  ^* 


(3.27) 


respectively.  The  Neumann  boundary  condition  for  the  nitride  interface  remains  un¬ 
changed: 

^  =  0.  (3.28) 


Oxidant  difftision  is  a  mixed  boimdary  value  problem  involving  three  different  types 
of  boimdary  conditions.  The  source  distribution  technique  is  flexible  in  dealing  with 
such  a  problem.  We  observe  that  both  Eqs.  3.7  and  3.11  share  the  same  source  p,  thus 
either  one  of  them  may  be  used  on  different  parts  of  the  boundary  to  describe  the 
problem,  depending  on  whether  the  boundary  conditions  are  Dirichlet  or  Neumann. 
A  linear  combination  of  both  produces  the  Robin  boundary  condition. 
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3.2  Biharmonic  Problems 


The  oxide  motion  process,  be  it  elastic  deformation  or  viscous  flow,  can  be  shown 
to  be  a  boundary  value  problem  related  to  the  Laplace  problem.  The  solution  to  aji 
elastic  deformation  or  incompressible  viscous  flow  in  two-dimension  is  biharmonic  or 
Bi-Laplacian: 

=  0  (3.29) 

w'here  the  bi-Laplace  operator  is  given  as 


V*  =  (V*)* 

d*  d*  d* 

dx\  dx\dx\  dx\ 

for  two  dimensions. 


Like  a  Laplace  problem,  a  biharmonic  problem  is  also  a  BVP  whose  solution  is 
uniquely  determined  by  the  boundary  conditions.  The  specification  of  the  solution 
requires  a  pair  of  boxmdary  conditions  from  the  list;  \1?,  5(V*'I')/3n. 

Once  again,  one  must  be  careful  in  picking  the  proper  pair.  If  ^  is  not  uniquely 
specified  anywhere  on  the  boxmdary,  the  system  becomes  ill-defined  because  ^  +  c, 
where  c  is  constant,  is  also  a  permissible  solution. 

Laplace’s  equation  can  be  considered  a  special  case  of  the  biharmonic  problem.  If 
a  function  /  is  harmonic: 

vV  =  0, 

it  always  satisfies  the  bi-Laplace  equation: 


VV  =  V*(VV) 
=  V’O 
=  0 


A  Laplace  problem  may  also  be  solved  as  a  bihaxmonic  problem  by  using  =:  0  as 
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one  of  the  boTindary  conditions  on  T.  Since  V*'i'  is  Laplacian,  is  zero  everywhere 
inside  Cl,  thxis  'i'  also  satisfies  Laplace’s  equation. 

At  one  time,  nonlinear  oxidant  diffusion  was  implemented  as  a  biharmonic  prob¬ 
lem.  The  idea  is  that  nonhomogeneity  in  the  diffusion  coefficient  results  in  pseudo 
domain-sources.  By  assuming  the  pseudo  domain-source  density  to  be  Laplacizin,  the 
system  becomes  biharmonic,  hence  domain  calculations  can  be  avoided.  However,  only 
passable  results  were  obtained. 

3.2.1  Biharmonic  Kernel  Functions 


Biharmonic  kernel  functions  can  be  readily  generated  from  the  harmonic  counterparts. 
We  note  that  if  f{x,y)  is  harmonic,  then 


f{x,  y)  itself  is  also  biharmonic,  as  noted  earlier. 

B  f 

is  biharmonic,  since  V^(z/(x,y))  =  =  0. 

Similarly,  y/(x,  y)  is  biharmonic. 

B  f 

(x^+y*)/(i,y)  is  also  biharmonic,  because  V^(z*+y*)/(x,y)  =  V*[4/ +4(i^ 


+ 


There  zu’e  more  varieties  than  the  Laplace  system.  As  it  will  be  shown  later,  we  will 
only  dealing  with  this  pair: 


<f>\  =  zlogjrl-x 

=  x<^  —  X 

<i>B  =  yiogjrl 
=  y<f>. 
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where  the  term  x  is  added  to  <f>\  to  make  it  compatible  with  the  period  version.  The 
corresponding  periodic  version  is  given  by: 

00 

H  (x  +  2nir)log|r  +  2n7rxiI 

n=— 00 
oo 

yIogIr  +  2n7r^i| 

n=— 00 
=  !/^*- 


There  is  no  simple  functional  form  for  But  as  discussed  in  a  previous  article  [18], 
its  derivatives  do: 


ilog[2(cosh(y)  -cos(x))l  - 

<p-y^ 

dy 

1  sin(x) 

«y - ^ - 


(3.30) 


2  cosh(y)  —  cos(x) 
'dx 


(3.31) 


We  note  that,  for  both  periodic  and  aperiodic  versions,  the  first  and  second  deriva¬ 
tives  of  and  <f>  share  many  common  terms.  Thus  it  is  appropriate  at  this  point 

to  introduce  a  set  of  basis  functions  consisting  of 

^,2>  ®2^,2>  X2<f>,ii,  ^2^,12 


where  indicial  notation  is  used  and  the  asterisk  (*)  is  dropped  for  simplicity.  The 
subscript  denotes  derivatives  in  the  i  direction;  thus. 


^2^,12  =  y 


ay 

dxdy’ 


These  expressions  are  combined  in  various  ways  to  generate  the  desired  Green’s  func¬ 
tions  for  diffusion,  viscous  flow,  elastic  deformation  and  viscoelastic  flow. 
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3.2.2  Generic  BEM  for  Oxide  Motion 

We  are  working  with  three  oxide  motion  models.  Except  for  the  kernels,  their  bihar¬ 
monic  BEM  formulation  are  similar  in  most  respects.  For  this  reason,  we  will  present  a 
template  BEM  framework  here,  delaying  the  derivation  of  the  kernels  until  we  discuss 
the  mathematics  of  the  appropriate  motion  model. 

Biharmonic  BEM  resembles  Laplace  BEM  in  many  ways.  The  displacement  vector 
u  parallels  the  scale  potential  The  s-urface  traction  p  is  similar  to  d$/dn.  The 
stress  tensor  Cij  is  analogoxis  to  V$.  There  are  also  the  direct  and  indirect  methods 
for  biharmonic  BEM.  Somigliana’s  formula  parallels  Green’s  formula.  The  vector 
potential  method  described  below  is  an  extension  to  the  scalar  source  used  in  the 
oxidant  diffusion  modeling. 

Vector  Potential  Method 

In  indicial  notation,  the  displacement  field  produced  by  boundary  sources  is  given  by 

u<(p)  =  ^  p,(q)«?,(q  -  p)dr  (3.32) 

where  ct’s  axe  the  “loading”  sources  and  u-y  is  the  displau:ement  in  the  x,-  direction  due 
to  a  loading  source  in  the  Xy  direction. 

Similarly,  the  formula  for  stress  is 

(p)  =  ^  pt(q)<^oi(q  -  p)dr.  (3.33) 

Using  the  definition  p,-  =  nya.y  (Eq.  2.16),  we  obtained  the  statement  for  surface 
traction: 

Pi(p)  =  ^  Pi(q)p<,(q  -  p)dr.  (3.34) 

Conceptually,  py  defines  loading  forces  at  the  boundziry,  as  shown  in  Fig.  3,2. 
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Figure  3.2:  Vector  Potential  Formulation.  The  source  density  p,-  defines  loading  forces 
applied  at  the  boundary  F. 

Somigliana’s  Formula 

Soznigliana’s  formula  is  the  foundation  for  the  direct  method  on  elasticity.  It  is  of  the 
form: 

«.(p)  =  ^  <;p,(p  -  q)  +  p*,«;(p  -  q)dr  +  ^  u;y(p  -  q)6,{p  -  q)dn  (3.35) 

where  the  normal  direction  n  (embedded  in  P”)  refers  to  the  dummy  boundary  point 
q.  As  we  can  see,  the  procedure  to  obtain  stress  information  for  interior  points  is 
rather  cumbersome  -  it  requires  the  derivatives  of  u,-,  and  that  means  a  new  set  of 
kernel  functions:  u.-y*  auid  p,*y*. 
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3.2.3  Incompressible  Slow  Viscous  Flow 


In  solving  an  incompressible  flow  problem,  it  is  customarily  to  express  the  velocity 
vector  in  terms  of  a  scale  parameter  called  stream  function.  The  x  and  y  components 
of  the  velocity  vector  v  =  v^x  +  v^y  is  related  to  the  stream  function  ^  by 


Vx 

v„ 


dy 

dx 


(3.36) 

(3.37) 


By  this  arrangement,  v  always  satisfles  the  incompressibility  equation: 


V-v 


d^_d^ 
dxdy  dxdy 
0, 


provided  ^  is  a  smooth  function  in  the  domain  of  interest. 


It  is  known  that  the  stream  function  for  incompressible  viscous  flow  is  biharmonic 


[43]: 


and  that  the  hydrostatic  pressure  is  a  harmonics  conjugate  function  to  the  vorticity 
u  = 


d(j}  dP 

^  dx  dy 

_  _£P 
^  dy  dx' 

Thus  the  stress  tensor  is  given  by: 


(3.38) 


(3.39) 


Now  we  select  —<pg  and  <f>\  as  Green’s  functions  corresponding  to  force  loading  in  the 
xi  and  3:2  directions  respectively.  The  complete  set  of  kernel  functions  are  given  below: 


Vii  =  —{X2<f>fi  +  <f>) 


(3.40) 
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Viz  =  (3.41) 

Cm  =  — 2(a;2^,ij  +  ^<f>,i)r]  (3.42) 

Cii2  =  2{x2(f>,ii  -  (3.43) 

<ti2j  =  2x2<?5.i2»7  (3.44) 

V21  =  X2^.i  (3.45) 

Vzz  —  {xz4>,z  —  <t>)  (3.46) 

<7211  =  2(X2^,11  -  <^,2)»;  (3.47) 

<7212  =  2x2<f>.izV  (3.48) 

<7222  =  ~2(X2<^,11  4*  <^,2)^7  (3.49) 


where  the  first  subscript  digit  denotes  the  direction  of  the  loading  force. 

3.2.4  Elastostatic  Deformation 

According  to  Papkovich  and  Neuber  [38],  a  general  solution  u  that  satisfies  the  Cauchy- 
Navier  equation  (Eq.  2.32)  can  be  generated  from 

u(p)  =  h(p)  -  (p  .  h(p)  +  /(p))  (3.50) 

where 

V’h(p)  =  0,  (3.51) 

V’/{p)  =  0.  (3.52) 

If  we  let  /  =  0  and  h  =  lp|“^xi,  we  obtain  the  fimdamental  solution  for  the 
elastostatic  problem.  Also  known  as  Kelvin’s  solution,  it  consists  of  the  stress  and 
displacement  fields  generated  by  a  point  force  in  an  infinite  3-dimensional  medium. 
The  body  force  is  applied  at  the  origin  r  =  0  in  the  xi  direction,  as  it  is  evident  from 
the  fact  that  the  Laplacian  is  violated  by  a  Dirac  delta  function  centered  at  that  point: 

V*1p|"^Xi  =  — 47r^(p)xi.  (3.53) 
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To  get  the  solutions  for  loading  forces  applied  in  the  other  two  principal  directions,  we 
change  h  to  and  Ipj'^xs  accordingly. 

3.2.5  Elastostatic  Kernels 

The  two-dimensional  Kelvin’s  solution  is  obtained  by  setting  /  =  0  and  h  =  log(l^)xi 
or  log(l^)x2.  It  is  used  as  the  kernel  or  Green’s  function  for  elastostatic  BEM  appli¬ 
cations  [42,40,18].  In  indicial  notations,  the  displacement  and  stress  fields  generated 
by  point  forces  in  the  xi  and  Xj  directions  are  given  below: 


“XI 

= 

(3.54) 

“12 

= 

4(1 

(3.55) 

<^111 

= 

(3.56) 

<^112 

_ 

(3.57) 

a’i32 

= 

(3.58) 

“21 

= 

4(1 

(3.59) 

“22 

= 

4(1 

(3.60) 

<^211 

= 

(3.61) 

<7212 

= 

(3.62) 

<7222 

= 

4(1 -■,)(! +^)  I"'" 

(3.63) 

Again,  the  first  digit  in  the  subscript  denotes  the  direction  of  the  loading  force,  as  in 
the  viscous  incompressible  flow  kernels. 


CHAPTER  3.  BOUNDARY  ELEMENT  FORMULATIONS 


66 


3.3  Laplace  Transform  and  Viscoelastic  BEM 

Elastostatic  deformation  is  a  static  problem  in  which  the  solution  does  not  change  if 
the  inputs  are  held  constant.  Incompressible  viscous  creeping  flow  is  a  steady-state 
problem  that  also  exhibits  similar  behaviors.  The  main  difference  between  them  is  the 
primary  variable  (that  may  be  used  as  an  input  condition)  is  displacement  for  elasticity 
and  rate  of  displacement  (i.e.  velocity)  for  flow.  On  the  other  hand,  viscoelastic  flow  is 
a  transient  problem.  Whether  constant  displacement  or  constant  velocity  is  specifled 
as  a  boundary  condition,  the  stress  field  still  changes  with  time. 

At  first  glance,  it  seems  that  linear  viscoelasticity  may  be  not  be  solved  in  a  reduced- 
dimension  manner  because  of  its  transient  natme.  However,  according  to  the  corre¬ 
spondence  principle  on  viscoelasticity,  a  viscoelastic  problem  may  be  modeled  using 
an  equivalent  elastostatic  problem.  In  simple  terms,  it  says  that  the  time-dependent 
stress-strain  behavior  in  the  bulk  is  related  to  the  boundary  conditions  in  a  simple 
manner. 

3.3.1  Laplace  Transform 

In  initial  value  problems  (IVPs),  the  solutions  are  defined  by  the  initial  conditions 
imposed  on  the  systems.  If  the  problems  are  given  in  terms  of  linear  differentia’ 
equations  with  constant  coefficients,  they  axe  known  as  linear  time-invariant  (LTI) 
systems.  The  Laplace  transform  method  is  a  convenient  and  indispensable  tool  for 
working  with  LTI  systems.  This  technique  converts  a  differentia  equation  into  a 
polynomial  function  of  a  transform  parameter.  After  the  unknown  coefficients  are 
calculated  algebraically,  the  inverse  transform  produces  the  solution  for  the  differential 
equation.  The  importance  of  the  Laplace  transform  is  highlighted  by  its  conversion  of 
time-domain  convolution  to  transform  space  multiplication.  This  simplication  assists 
the  analysis  of  complex  interconnected  subsystems. 
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The  definition  of  Laplace  transform  and  its  inverse  are  given  in  the  following  no¬ 
tation: 


Fis)  =  £{/(t)} 

=  /  fiT)e~*^dT 

(3.64) 

fit)  =  £-^{Fis)} 

1  /a-Hoo 

where  i  is  the  imaginary  number.  The  lower  limit  on  the  Laplace  transform  is  set  to 
0~  so  that  impulse  behaviors  at  f  =  0  in  f(t)  are  included  in  the  integration. 

Important  characteristics  relevant  to  our  application  are  given  below: 


£ia  +  b}  =  j:{a}  + £{5} 
£{da/dt}  =  s£{a}  +  il{a(0)} 
flit)  *  hit)  =  F^{s)  X  F^{s) 
/(O’*")  =  lim,_oo5F(s) 

/(co)  =  lim,-.os.F(s) _ 


Lineax 

Time  Differentiation 
Time  Convolution 
Initial  Value  Theorem 
Final  Value  Theorem 


where  the  definition  of  convolution  integral  is 

flit)  *  /2(0  =  /  fiir)f2it  -  r)dr. 

OO 


Actual  evaluations  of  trainsform  integrals  are  seldom  required  since  table  lookup  is 
more  convenient.  Given  in  Table  3.1  axe  a  list  of  formulae  useful  for  our  applications. 


3.3.2  Viscoelastic  Kernels 

As  stated  in  Chapter  2,  the  constitutive  equation  relating  stress  and  strain  is  usually 
given  in  time-derivative  form,  as  in  Eq.  2.42.  Using  the  Laplace  transform  technique, 
it  can  be  converted  to  a  polynomial  expression  of  the  transform  parameter  s: 

Pis)a  =  (3.65) 
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Inverse  Transform  F{s) 

Function  f(t) 

1 

-m 

1/s 

tt(t) 

lls^ 

u{t)t 

l/(s  +  a) 

e-‘“u(t) 

Table  3.1:  Laplace  Transform  Table 


where 

P  = 

t=0 

«=0 

Note  that  we  have  made  an  implicit  a;5sumption  that  the  system  is  initially  at  rest, 
ue.  c{t  =  0)  =  €{t  =  0)  =  cr(t)  =  €(t) . .  •  =  0. 

The  strategy  for  finding  the  viscoelastic  kernels  is  to  employ  expressions  already 
derived  for  linear  elasticity.  Since  most  of  them  are  stated  in  terms  of  Young’s  modulus, 
E,  and  Poisson’s  ratio,  u,  it  is  desirable  to  relate  the  Laplace  polynomials  to  these  two 
parameters.  The  relevant  constitutive  equations  given  in  Chapter  2  zu-e  reproduced 
here  for  comparison;  they  axe 

Sij  —  2G€ij 

Cm  — 

for  elasticity  and 

Pl^ij  — 

P2<^M  = 


for  viscoelasticity.  From  them,  the  following  relationships  are  obtained: 
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3K 


Since  G  =  E/2{l  +  i/)  and  K  =  E/3{l  —  2i/)t  the  expressions  for  E  and  1/  are  therefore 


obtained[36|; 


i;(s) 

I/(s) 


3<3i(s)^2(s) 

^2i(>s)P2(s)  +  2Pi(s)  <22(5) 

Pl(^)g2(^)-gl(3);>2(3) 

i2l(3);>2(5)+2Pi(s)fi2(5)’ 


(3.66) 

(3.67) 


With  these  modified  definitions  E(s)  and  i^{s),  the  solution  for  an  elastic  problem  can 


be  converted  to  a  viscoelastic  form  by  applying  the  inverse  Laplace  transform. 


This  conversion  procedure  also  applies  to  the  elastic  kernels  Eqs.  3.54-3.63  to  gen¬ 
erate  a  set  of  viscoelastic  kernels.  Notice  that  we  can  premultiply  all  these  expressions 
with  a  scale  factor  W  (s)  before  the  inverse  Laplace  transform.  This  effectively  modifies 
the  behavior  of  the  loading  force.  Hence  we  can  tailor  W  (s)  to  generate  a  solution  with 
desirable  characteristics  for  our  modeling  effort.  Typically,  we  set  either  the  displace¬ 
ment  or  the  force  (but  not  both)  at  the  load  point  to  follow  a  step  change  or  linear 
change  in  time.  The  drivmg  force  of  oxidation  is  the  oxide  growth  at  the  oxide-silicon 
interface.  As  a  first  order  approximation,  we  assiime  that  the  growth  rate  is  constant 
within  each  time  step,  therefore  the  injection  velocity  is  also  constant.  This  suggests 
that  a  constant-velocity  formulation  is  best  suited  for  our  application.  For  some  other 
applications,  a  constant-force  loading  may  be  more  appropriate. 


To  obtain  a  constant- velocity  loading,  we  note  that  at  the  load  point  (r -+  0),  Um 
(or  U222)  is  dominated  by  the  term  in  log  |fj,  which  has  a  coefficient  of 

4t/(s)  —  3 
4(l-i/(s))’ 

The  premultiplier  is  set  to 


4(1  -  i/(s))  1 


(3.68) 


41/(5)  —  3  5* 

to  convert  the  coefficient  to  1/s*,  which  is  the  transform  of  a  constant  velocity  {du/dt  — 
c).  Putting  it  in  Eqs.  3.54-3.63  and  taking  the  inverse  Laplace  transform  results  in 
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the  following  expressions: 

ttii  =  7Kc(l  -  -  (x2^,2  +  (l>)t  (3.69) 

tti2  =  “7iifa(l  “  2i/)x2^,1  +  X2^,lt  (3.70) 

am  =  [-Kai{2x2<f>,i2  +  Z(f>^i)  —  Kfi{x2(f>,i2  +  5<f>,i)]E  (3-71) 

<^112  =  [Ka{2X24>,ll  —  ^,2)  +  Kfi{X24>,ll  —  4(i>,2)]E  (3.72) 

ai22  =  [■K^a(2X2<^,12  “  ^,l)  + +  3^,l)]i?  (3.73) 

U21  =  -7iifa(l  —  2i/)X2<^.l  +  X2^.lt  (3-74) 

U22  =  —7Ka{l  —  2i/)x2^,2  +  (®2^,2  “  4)^  (3.75) 

0'211  =  (■^^a(2X2^,U  “  3^,2)  +  ■K^s(a^2<^.U  +  2^,2)].^  (3.76) 

C212  —  [-^Ca  (2X2<^, 12  +  ^.1)  +  ■^i9(^2<^, 12  —  3^, l)]i^  (3.77) 

a222  =  (•^^<*(222^,11  +  ^.2)  +  ■^^/9(3:2^,n  +  ^^,2)]-^  (3.78) 


where 

=  ffll  - 

3(3  —  4u)rj 

^ - 

2(1  +  v)‘n 
=  ~E  ~ 

_  ]l 

~  G' 

Ka  and  are  analogous  to  the  charging  of  the  capacitor  C  in  a  simple  R-C  circuit. 
Like  the  capacitor  that  dischsxges  when  the  applied  voltage  is  removed,  and 
change  to  decaying  exponential  functions  with  time  constants  r*  and  respectively 
when  the  loading  is  stopped. 


In  general,  the  two  relaxation  time  constants  zure  different.  The  second  one,  r^,  is 
readily  identified  as  shear  relaxation.  The  origin  of  Ta  is  less  clear;  probably  it  has  to 
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do  with  combined  volume-sheax  relaxation  as  evident  from  the  fact  that  it  becomes 
the  same  as  the  shear  relaxation  time  when  the  material  is  incompressible  {u  =  0.5). 
Stress  fields  associated  with  and  Kf  satisfy  the  force  conservation  equation  Fq.  2.17 
independently.  This  condition  arises  from  the  fact  that  Ka  and  decay  differently. 

We  expect  that,  as  the  ratio  rj/G  approaches  oo  or  0,  the  viscoelastic  kernels  reduce 
to  that  of  elastic  deformation  or  viscous  incompressible  fiow  respectively. 

For  the  first  case,  the  limit  t]/G  -*  oo  vs  obtained  by  keeping  G  at  a  finite  value 
and  letting  rj  go  to  infinity.  The  exponential  functions  are  approximated  with  Taylor 
series  and  the  following  expressions  are  obtained 

=  —!—t 
3(3  -  4u) 

E^  ^  E  ^ 

"  2(1 +  1/)*' 

for  the  appropriate  terms  in  Eqs.  3.69-3.78.  In  this  regime,  the  viscoelastic  kernels  are 
the  elastic  kernels  given  in  Eqs.  3.54-3.63  scaled  by  4(1  —  i/)tf{4i/  —  3). 

In  the  second  case,  the  limit  ijjG  -*  0  is  reached  by  fixing  77  at  a  finite  value  and 
letting  G  go  to  infinity.  For  any  nonzero  t,  no  matter  how  small,  the  exponentials 
are  zero.  We  find  that  the  viscoelastic  kernels  simpliiy  to  the  viscous  incompressible 
fiow  kernels  (Eqs.  3.40-  3.49).  Note  that  1/  does  not  appesir  in  the  final  expressions. 
The  rezison  is  that  the  elastic  mod;ili  are  set  to  be  infinitely  large  in  this  simplication 
(G  -+  00) ,  th\is  the  material  is  incompressible  irrespective  of  u. 

The  viscoelastic  kernels  are  different  from  those  of  elasticity  and  viscou  fiow  in  that 
they  have  an  additional  parameter  t.  This  pau'ameter  is  set  to  the  time-step  size  At  in 
numerical  simulations. 


As  we  can  see,  the  constant-velocity  (linear-displacement)  loading  viscoelastic  ker- 
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nels  approach  is  a  flexible  formulation.  Because  the  stress  relaxation  time  can  be 
pushed  to  the  extremes,  this  viscoelastic  formulation  is  also  suitable  for  modeling 
elastostatic  deformation  and  viscous  incompressible  flow. 


3.4  Modeling  Nonlinear  Bulk  Behaviors 


The  integral  equation  in  the  BEM  is  nothing  more  than  a  continuous  summation  of 
contributions  from  boundary  sources.  For  this  superposition  technique  to  be  valid,  the 
domain  has  to  be  homogeneous,  t.e.  bulk  parameters  such  as  diSusivity  and  viscosity 
have  to  be  uniform  throughout  the  bulk,  as  originally  assumed  in  the  derivations  of 
those  Green’s  functions.  It  has  been  mentioned  earlier  in  Chapter  2  that  stress  alters 
those  parameters,  rendering  the  domain  nonhomogeneotis.  This  section  describes  a 
formalism  for  handling  the  nonlinear  effects  with  the  BEM. 


The  concept  comes  &om  the  observation  that  both  the  diffusion  flux  and  stress- 
strain  equations  are  of  the  form 


A  =  cB. 


(3.79) 


where 


Diffusion 

Elasticity 

A 

F 

a 

B 

VC 

e 

c 

D 

E 

For  notational  simplicity,  the  strain  and  stress  parameters  are  assumed  to  be  scalars. 

The  similarity  suggests  that  both  processes  can  be  treated  in  a  unified  manner.  For 
generality,  we  will  discuss  the  issues  in  terms  of  the  representative  pareuneters  A,  B, 
and  c. 

Nonhomogeneity  or  nonlinearity  implies  the  parameter  c  is  not  uniform  within  the 
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domain,  c  can  be  rewritten  as 

c  =  coA:(x,.4,-*-)-  (3.80) 

As  shown,  the  scale  factor  k  may  be  a  function  of  many  parameters,  including  A,  B 
and  position  x.  To  conform  with  the  linear  system  format,  Eq  3.79  is  rewritten  as 

A  =  CoB  +  A°,  (3.81) 

where  A°  is  co(fc  —  l)B.  In  a  normal  setting,  A°  is  interpreted  as  a  initial  condition 

of  A  that  does  not  have  a  associated  B^  value.  If  B  is  expressed  in  terms  of  A,  this 

alternate  statement  is  obtained: 

B  =  -^A  (3.82) 

=  —A  +  B°  (3.83) 

Co 

where  B°  is  the  initial  condition  of  B. 

One  can  find  familiar  situations  in  mechanics  where  there  is  A  without  the  accom¬ 
pany  B,  and  vice  versa.  For  instance,  consider  a  freely-supported  material  undergone 
thermal  expansion:  it  sees  a  nonzero  thermal  strain  But  it  has  zero  thermal  stress 
becatise  it  is  freely  supported.  On  the  other  hand,  if  the  material  is  clamped  on 
all  sides,  then  its  ^  is  zero  because  it  is  not  free  to  expand.  Consequently,  its  is 
non-zero.  It  is  said  to  have  initial  stress.  For  the  difiusion  problem,  similar  analogies 
2ire  hardy  to  come  up  with.  We  note  that  in  insulating  materials,  the  increase  in 
dielectric  constant  (beyond  cq)  is  due  to  electronic  polarization.  If  we  can  freeze  the 
polarization  when  we  remove  the  applied  field,  we  will  see  a  residual  built-in  field.  Be¬ 
cause  dielectric  constant  in  potential  is  analogous  to  diffusion  coefficient  in  diffusion, 
we  coin  the  term  “built-in  field”  as  the  initial  condition  for  Laplace  problems. 

3.4.1  Viscoplastic  BEM 

This  unified  framework  for  nonlinear  diffusion  and  viscoelastic  flow  evolves  from  consid¬ 
erations  on  elastoplastic  BEM  method.  Elastoplasticity  is  the  permanent  deformation 
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a 


OR  e 


Figure  3.3:  Elastoplasticity.  Beyond  the  point  P,  the  system  is  nonlineax.  The  effective 
E'  is  different  from  E  of  the  lineax  regime. 

in  an  elastic  material  due  to  large  stress  or  strain.  Typically,  such  an  effect  b  observed 
only  when  stress  exceeds  the  yield  limit  of  the  material.  Fig.  2.5  from  Chapter  2  b 
reproduced  here  as  Fig.  3.3  to  show  a  representative  plastic  stress-strain  curve 

An  unmodified  BEM  can  only  model  the  initial  linezir  regime,  as  defined  by  OP 
in  the  figure.  Beyond  the  point  P,  the  stress-strain  relationship  b  nonline^.  At  Q, 
the  effective  E'  b  given  by  cq/cq.  Because  stress  (or  strain)  b  not  likely  to  uniform 
throughout  the  domain,  different  regions  will  have  different  effective  E'.  Thus  the 
domain  b  not  homogeneous. 

The  three  common  techniques  for  nonlinear  elasticity  are  the  modified  body-force, 
initial  stress,  and  initial  strain  formulations.  Only  the  concepts  of  initial  strain  and 
initial  stress  will  be  described  here  as  they  relevant  to  our  work.  For  details  on  all 
these  methods,  readers  are  referred  to  the  book  authored  by  Banerjee  and  Butterfield 
3  In  this  discnsaion,  we  are  not  concerned  with  how  the  system  returns  to  the  unstressed  state. 
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[40].  The  treatment  by  Brebbia  et  al  [42]  Is  not  as  in-depth  as  the  former. 


In  the  “initial  stress”  formulation,  the  nonlinear  stress  components  is  expressed  as 


[40]: 


(3.84) 


where  a°  is  the  initial  stress,  Ct  is  the  elastic  stress,  and  a**'  is  the  elastoplastic  stress 
that  we  seek.  The  equilibrium  equation  then  becomes 


dxj  dxj ' 

As  far  as  the  elastic  component  is  concerned,  there  is  a  pseudo  body  force  due  to  initial 
stress  a°.  The  displacement  formula  is  modified  to  include  this  contribution: 

t^(p)  =  ^  py<y(p  -  q)dr  +  ^  -  q)<in.  (s.se) 

The  statement  for  elastoplastic  stress  is  accordingly  given  by 

“*  (3-87) 

Note  that  the  stress  expression  has  a  term  in  that  has  no  counterpart  in  the  displace¬ 
ment  formula.  The  domain  integrals  in  both  expressions  are  usueilly  converted,  using 
the  divergence  theorem,  to  a  different  form  that  utilizes  a?-  as  the  domain  source, 
instead  of  5or?  /5xy.  With  this  transformation,  we  avoid  taking  derivatives  of  5cr?y, 
which  may  have  to  be  done  numerically.  However  there  is  no  need  to  do  so  in  our 
particular  formulation;  hence  we  stay  with  the  original  expressions. 


In  the  initial  strain  formulation,  we  modify  strain  instead  of  stress: 

^ij 


(3.88) 


where  is  the  initial  strain,  e,y  is  the  total  strain,  and  €*y  is  the  elastic  strain.  Stress 
is:  given  in  terms  of  the  elastic  strain: 
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(b) 

Figure  3.4:  Modeling  nonlinear  elasticity  with  (a)  the  initial  stress  approach  and  (b) 
the  initial  strain  approach. 
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The  initial  strain  6°y  is  then  used  in  domain  integrals  similar  to  those  for  £r,°. 

Fig.  3.4  illustrates  how  nonlinear  elasticity  is  modelled  with  the  initial  stress  and 
initial  strain  formulations.  To  reach  the  point  Q  defined  in  Fig.  3.3,  the  initial  stress 
approach  uses  a  vertical  translation  (stress  without  strain  -  initial  stress)  in  conjunc¬ 
tion  with  a  linear  deformation.  Mathematically  this  is  given  by  Eq.  3.81.  Likewise, 
the  initial  strain  formulation  employs  a  horizontal  translation  (strain  without  stress  - 
initial  strain)  and  a  linear  deformation,  as  represented  by  Eq.  3.83. 

3.5  Nonlinear  Diffusion  BEM 

The  nonlinear  diffusion  follows  the  initial  stress  method  closely.  The  oxidant  flux 

F  =  -DVC 
=  — V  C 

=  F^  +  F°  (3.89) 

is  split  into  a  linear  (or  homogeneous)  part 

=  -PoVC 

and  a  built-in  field 

F°  =  Dq{1  -  ki)VC 
where  kd  is  the  diffusivity  scale  factor. 

The  fl\ix  conservation  is  on  F: 

V.F  =  V.F^  + V-F° 

Hence  the  linear  part  sees  a  pseudo  domain  source  due  to  the  built-in  field: 

p  =  -V-F^ 

=  V-F°. 


(3.90) 

(3.91) 

(3.92) 
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The  integral  equation  for  C  now  becomes 

C{p)  =  /  ^(q)^-(q-p)dr+  /(V,-r°(q))^-(q-p)dn  (3.93) 

Jr  Jfi 

where  the  subscript  in  V,  indicates  the  operation  is  done  with  respect  to  the  dummy- 
parameter  q,  rather  than  with  p.  The  integral  equation  for  flux  is  accordingly  given 
by 

F  =  DoVC  +  F° 

=  Do  /■  (.(q)V,i-Cq  -  p)dr  +  Do  /  (V,  ■  F“(q))V^'(q  -  p)dn  +  F” 

Jr  J(i 

=  Dokj  [1^  /.(q) V*-(q  -  p)dr  +  fjV,  ■  F°(q))V«-(q  -  p)dn]  (3.94) 


1 


3.6  Nonlinear  Viscoelastic  BEM 


Our  nonlinear  viscoelastic  BEM  formulation  is  different  from  the  elastoplastic  or  vis¬ 
coplastic  formulations  because  of  its  use  of  the  Laplace  transform  technique.  In  the 
conventional  approach,  the  deviatoric  components  in  stress  (or  strain)  is  singled  out 
for  plastic  deformation  treatment,  leaving  the  spherical  components  unchanged.  In  our 
problem,  we  are  only  concerned  with  the  effect  of  stress  on  the  viscosity.  The  viscosity 
parameter  appears  in  the  exponential  functions  Kx  Kb  of  tihe  viscoelastic  kernels 
Eqs.  3.69 —  3.78.  Given  that 

»7  =  (3.95) 


where 


(^sYrjilkT 


kr,  =  exp  j  .  (3.96) 

^  \  kT  J  smh.{asY,)2fkT)  ' 

as  described  in  Chapter  2,  we  can  expressively  show  the  influence  of  kr,  by  writing  the 
kernels  of  displacement  and  stress.  For  simplicity,  we  introduce  two  sets  of  variables 
to  represent  the  effective  scaling  of  k,,  on  the  displacement  and  stress: 

tt,(p)  =  ^Pi(q)tt*,.(q-p,At,A,)dr 
=  ^u.(p)ttf(p) 
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where  the  superscript  L  denotes  the  solution  to  the  linear  system  in  which  k„  =  0. 


The  nonlinear  stress  must  satisfies  the  force  equilibrium  condition: 

dXj 

-  ’axj' 


This  results  in  a  pseudo  body  force  of 


6(  =  aafjdxj 

= 

kg  **  dxj 

for  the  linear  part. 

The  presentation  so  far  has  been  analytical  and  general.  How  well  this  formulation 
performs  depends  critically  on  the  numerical  approximation  technique. 


(3.97) 


3.7  Incompressible  Materials 


In  certain  numerical  tests  on  ring  structmes,  the  loading  force  is  arranged  in  a  radial 
direction,  as  shown  in  Fig.  3.5.  An  anomalous  behavior  was  found  when  Poisson’s 
ratio  is  set  to  0.5,  *.c.  when  the  material  is  incompressible.  Outside  the  ring,  no 
displacement  and  stress  is  seen  (except  for  residual  numerical  errors).  Inside  the 
ring,  only  a  high  pressure  is  observed,  the  displacement  remains  zero.  The  inability 
to  generate  a  radial  symmetric  solution  reveals  that  Kelvin’s  solution  that  we  have 
been  using  is  deficient  —  it  cannot  model  a  certain  legitimate  behavior  or  geometric 
configxiration. 
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Figure  3.5:  Radial  loading  forces.  Neither  stress  nor  displacement  is  observed  outside 
the  ring. 

This  is  problematic  for  modeling  any  annulus  structure  because  it  has  an  outer 
circular  boundary  and  an  inner  circular  boundary.  The  inability  of  the  inner  boundary 
sources  to  influence  the  domain  behavior  means  the  matrix  solution  is  singular  or  nezurly 
so,  producing  very  poor  results.  Granted,  one  would  use  a  more  realistic  Poisson’s  ratio 
of  0.125,  for  instance.  However,  as  we  noted  in  a  previous  section,  the  material  will 
automatically  becomes  incompressible  as  the  viscosity  decreases. 

Note  that  the  inner  boundary  needs  not  be  circular  for  singular  behavior  to  show 
up  -  abnormal  results  are  also  obtained  when  the  inner  boxmdary  is  rectangular. 

The  physical  cause  of  the  singular  behavior  can  be  understood  by  examining 
Fig.  3.6.  Shown  in  part  (a)  is  an  incompressible  annulus  structure  that  is  to  be 
modeled.  The  inner  hole  zmd  the  exterior,  which  are  not  shaded,  are  of  a  different 
material  that  is  totally  compressible  -  free  space.  The  way  the  indirect  BEM  han- 
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Figure  3.6:  Source  of  incompressible  singularity,  (a)  Actual  problem,  (b)  In  BEM 
implementation,  the  hole  in  the  center  is  filled  with  the  same  incompressible  material. 

die  this  problem  is  illustrated  in  part  (b).  The  whole  space  is  filled  with  the  same 
incompressible  material,  as  required  by  Kelvin’s  solution.  The  BEM  applies  loading 
force  on  the  two  circular  lines  that  map  to  the  actual  boundaries  in  (a)  so  that  the 
"boundary  conditions”  on  these  lines  match  the  real  ones.  The  problem  now  is  that 
the  hole  is  filled  with  an  incompressible  material.  The  Kelvin’s  solution  only  produce 
distortion,  it  does  not  introduce  any  new  material  or  otherwise  change  the  volume. 
Hence  it  cannot  change  the  area  of  the  hole,  and  this  behavior  is  not  compatible  with 
the  real  situation.  Thus  a  special  solution  is  needed  in  order  to  be  able  to  alter  the 
size  of  the  hole. 

The  special  solution  is  an  elastic  potential  function: 

=  ^log(i*  +  y*). 


CHAPTER  3.  BOUNDARY  ELEMENT  FORMULATIONS 


82 


The  displacement  field  is  given  by 


and  hence  the  stress  is 


X 


'11 

_0 


X*  +  +  y2‘ 


=  2G 


—  x^ 

(PTy^ 


cl,  =  0 


'22 


=  2G- 


x»-y* 


(3.98) 

(3.99) 


(3.100) 

(3.101) 

(3.102) 


(x*  +y*)*‘ 

As  we  can  see,  there  is  no  spherical  deformation  or  stress  associated  with  this  function, 
except  at  the  origin: 


^11  +  ^22  =  V  •  Tl(p)  (3.103) 

=  22r^(p).  (3.104) 

where  there  is  a  source  of  material.  The  corresponding  viscoelastic  functions  are 
obtained  by  scaling  the  displacement  field  with  t  and  replacing  2G  in  the  stress  ex¬ 
pressions  with  2r][l  —  sxp(— Gt/f])]. 

To  include  this  special  feature  in  the  numerical  solution,  an  additional  condition  on 
the  system  is  needed.  After  experimenting  with  several  options,  it  is  decided  that  the 
best  constraint  is  to  specify  the  net  hydrostatic  pressure  generated  by  the  boundary 
sources  be  zero  at  the  origin.  With  this  enhancement,  it  is  found  that  the  solutions 
are  very  well  behave. 


Chapter  4 

Numerical  Solutions 


Solving  a  boundary  value  problem  with  the  BEM  entails  the  following  steps.  First, 
the  boundary  is  suitably  divided  into  segments  or  elements.  The  sources  are  assumed 
to  take  a  certain  shape  distributions  on  these  elements.  A  matrix  is  set  up  to  define 
the  sources  in  terms  of  the  known  boundary  conditions.  After  the  source  values  are 
obtained  from  the  matrix  solution,  the  unknown  boundary  and  interior  parameters 
can  then  be  calculated. 

In  modeling  thermal  oxidation,  the  solution  of  oxidant  diffusion  is  first  computed 
to  determine  the  oxide  growth  rate  along  the  silicon  interface.  The  injection  of  newly- 
created  oxide  into  the  bulk  is  calculated  from  the  growth  rate  and  desired  time-step 
size.  This  information  is  then  fed  into  the  fiow  matrix  equation  that  may  represent 
oxide  motion  as  elastic  deformation,  incompressible  viscous  flow,  or  viscoelastic  flow. 
From  the  flow  solution,  the  oxide  displacement  at  the  free  surface  is  deduced.  At  this 
point  the  motion  pattern  of  all  boundary  points  are  known,  hence  their  locations  can 
be  updated.  A  new  step  can  then  begin.  This  loop  continues  until  the  desired  oxide 
thickness  or  oxidation  time  is  reached. 

When  the  mechanical  interaction  of  oxide  with  the  nitride  layer  and  effects  of  stress 
are  considered,  the  numerical  procedures  become  more  complex  aind  intertwined  -  ma- 
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trbc  solutions  are  more  complicated,  domain  calculations  are  required,  and  iterations 
are  embedded  in  other  iterations.  The  flow  chart  described  in  the  previous  pzuragraph 
serves  to  give  a  clear  outline  of  the  underlying  activities.  The  presentation  of  this 
chapter  does  not  necessarily  follow  the  order  of  program  execution.  Rather,  the  de¬ 
scription  is  organized  aroimd  logical  groupings  of  functionalities.  The  segmentation  of 
boundary  and  domain  will  flrst  be  discussed.  After  issues  of  line  integration,  matrix 
synthesis  and  solution  will  be  dealt  with.  Methods  for  solving  nonlinear  equations  will 
also  be  presented. 


4.1  Boundary  and  Domain  Segmentation 


The  flrst  step  in  BEM  numerical  approximation  is  to  discretize  the  curved  boundaries 
by  subdividing  them  into  small  segments  or  elements.  This  subdivision  serves  two 
purposes.  The  piece-wise  straight-line  approximation  makes  integral  calculations  more 
tractable.  It  also  provides  a  suitable  mechanism  to  prescribe  a  shape  function  for  the 
source  distribution  on  each  segment. 

Three  basic  boundary  segmentation  techniques  are  shown  in  Fig.  4.1;  they  are 
(a)  constant,  (b)  linear,  and  (c)  quadratic.  In  the  constant  element  representation, 
the  collocation  point  at  which  the  boundary  condition  is  specified  is  located  in  the 
middle  of  an  segment.  For  the  indirect  method,  the  “constant”  term  refers  to  the 
fact  that  the  soxirce  density  along  a  segment  is  uniform.  Its  values  takes  a  step  jump 
from  one  segment  to  another,  as  depicted  in  Fig.  4.2a.  For  the  direct  method,  the 
“constant”  term  refers  to  the  weighting  function  on  the  boundeiry  condition.  In  the 
linear  element  configuration,  the  collocation  points  are  the  endpoints  of  the  segments. 
Source  density  is  assumed  to  vary  linearly  from  one  collocation  point  to  the  adjacent 
one,  as  shown  in  Fig.  4.2b.  Both  “constant”  and  “linear”  segments  are  straight  lines. 
In  the  quadratic  case,  a  boundary  element  is  a  curve  defined  by  3  nodes  points.  Since 
there  are  3  collocation  points,  it  is  therefore  possible  to  make  a  quadratic  fit  for  the 
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Figure  4.1:  Bo\indary  segmentation  techniques:  (a)  constant  element,  (b)  linear  ele¬ 
ment,  and  (c)  quadratic  element.  The  “nodes”  are  the  collocation  points. 

source  distribution. 

In  principle,  the  smoother  the  distribution  approximation  is,  the  better  the  results 
will  be.  However,  only  the  constant  element  is  practical  for  the  indirect  method; 
the  other  two  are  actually  only  relevant  to  the  direct  method.  The  reason  is  that 
these  segmentation  schemes,  with  the  exception  of  the  quadratic  approach,  are  just  a 
polygonization  of  a  geometry.  Where  two  boundary  segments  are  joined  together,  the 
surface  curvature  is  infinite.  In  other  words,  we  have  a  Liapunov  surface  that  is  not 
smooth  [38].  Whenever  the  surface  orientation  of  a  boundziry  point  is  discontinuous, 
the  indirect-method  solutions  are  singular.  No  meaningful  boundary  condition  can 
be  imposed  or  extracted  if  a  collocation  point  is  placed  there.  Therefore  collocation 
points  are  set  at  the  center  of  each  segments,  as  in  the  constant  element  scheme. 


A  linear-element  approximation  had  been  tried.  To  avoid  stirface  normal  discon¬ 
tinuity,  the  edge  joining  two  boimdary  segment  was  replaced  by  an  arc.  The  arc 
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Figure  4.2:  Source  distributions:  (a)  constant  element,  and  (b)  linear  element. 

curvature,  which  is  the  reciprocal  radius  of  a  circle,  is  of  a  predefined  value.  An  at¬ 
tempt  to  improve  the  “mid-point”  collocation  method  had  also  been  experimented. 
The  source  distribution  was  made  to  vary  linearly  from  one  collocation  point  to  the 
next.  Both  approaches  were  abandoned  because  the  minor  improvement  in  accuracy 
could  not  justify  the  increased  complexity  in  the  solution  technique.  This  point  is 
particularly  applicable  to  our  oxidation  modeling  effort  because  we  aire  trying  to  ex¬ 
plore  the  usability  of  the  boundary  element  technique,  rather  than  to  tune  it  to  get 
improved  computation  performance. 

4.1.1  Classes  of  geometries 

As  mentioned  in  Chapter  3,  it  is  desirable  to  mirror  a  simulation  structure  in  such  a 
way  that  the  upper  boundziry  (i.e.  the  free  surface  and  the  silicon  nitride  mask)  never 
touches  the  lower  boundary  (i.e.  the  silicon  substrate).  This  is  to  avoid  artificiztl 
boundary  closure  that  can  introduce  local  numerical  errors  in  regions  next  to  the 
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artificial  wall.  Two  types  of  mirroring  formats  are  supported.  First,  the  mirroring  of 
a  structure  at  the  left-  and  right-hand  edges  results  in  periodic  symmetry.  This  form 
requires  periodic  kernels.  The  other  mirrors  one  horizontal  edge  and  one  vertical  edge 
to  obtain  a  2-fold  symmetry;  it  uses  the  normal  non-periodic  kernels.  A  vajiation  on 
the  2-fold  symmetry  is  the  4-fold  one  that  has  an  additional  folding  at  a  45®  line.  All 
these  configurations  are  depicted  in  Fig.  4.3. 

Reduction  in  number  of  unknowns  is  achieved  by  exploiting  symmetry  features  of 
the  structures.  The  solution  is  only  calculated  for  half  a  period,  a  quadrant,  or  hailf  a 
quadrant  respectively.  Typically,  a  total  of  40  to  60  segments  zire  used  for  periodic  and 
2-fold  symmetric  structures,  and  20  to  30  for  the  4-fold  symmetric  case.  The  lower 
limit  is  due  to  accmacy  consideration  of  the  discrete  approximation  to  continuous 
quantities  while  the  upper  limit  is  constrained  by  the  computational  speed. 

Depending  on  the  nature  of  the  kernel  fimctions,  the  source  distributions  are  either 
symmetric  or  antisymmetric  with  respect  to  a  symmetry  line.  The  kernels  for  the 
oxide-flow  problem  work  in  pairs  to  generate  a  loading  force.  The  ones  specifying 
loading  forces  in  the  x  direction  produce  a  source  distribution  that  is  antisymmetric 
with  respect  to  a  vertical  fold  line  (o’(x)  =  — a(— z))  but  symmetric  with  respect  to  a 
horizontal  fold  line  (<r(y)  =  <t(— y)).  The  converse  is  true  for  the  y— loading  kernels. 
Moreover,  in  the  4-fold  symmetry  case,  the  source  distribution  of  a  “x"  kernel  is 
mirrored  to  that  of  a  “y”  kernel  along  the  45°  line.  For  the  oxidant  problem,  the 
kernel  is  scalar,  therefore  the  source  distribution  is  always  symmetric  with  respect  to 
the  fold  lines. 

The  periodic  symmetric  geometry  is  used  for  analyzing  LOCOS  and  SWAMI  struc¬ 
tures.  The  trench  structure  has  a  long  vertical  perimeter,  solving  it  as  a  periodic 

structure  will  require  many  boundary  segments  and  slow  down  the  computation  speed 

.'1^ 

significantly.  Therefore  we  take  advantage  of  the  vertical  wall  and  cut  it  in  the  middle 
into  an  upper-half  and  lower-half,  and  solve  each  of  them  as  a  2-fold  symmetric  object. 
The  4-fold  symmetry  is  primarily  used  for  modeling  Kao’s  experiments  on  cylindrical 
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structures. 

4.1.2  Interior- Cell  Generation 

It  is  possible  to  subdivide  a  nonhomogeneous  domain  into  piece-wise  homogeneous 
pzirtitions  and  solve  them  as  a  group  of  interconnected  systems.  This  subregion  tech¬ 
nique  is  effective  when  the  bulk  parameters  are  different  in  a  few  large  regions.  It  is 
used  for  modeling  the  mechanical  interaction  between  oxide  and  nitride  bodies.  How¬ 
ever,  it  is  not  suitable  for  modeling  stress  effects  because  many  partitions  may  be 
needed  to  model  the  variations  in  bulk  parameters  accurately;  this  complicates  the 
matrix  solution  techniques.  Thus  perturbation  techniques  based  on  interior  cells  are 
explored. 

Because  of  the  ever-changing  domain,  cell  partitions  are  not  specified  in  the  input 
data  but  are  determined  by  the  program  on  the  fly  at  every  time  step.  To  write  a 
fully  automated  generator  is  a  difficult  task;  only  a  rudimentary  scheme  has  been  used. 
This  method  takes  advantage  of  the  fact  that  upper  and  lower  boundziry  segments  axe 
aligned  and  specified  in  pairs.  Where  interior  cells  are  needed,  it  creates  a  coliimn 
joining  the  segment  pairs,  as  illustrated  in  Fig.  4.4a.  Then  it  looks  at  the  dimension 
of  the  column  to  determine  the  number  of  cells  that  should  be  created,  using  this  rule: 

effective  height 
effective  width  ’ 

where  the  effective  height  is  the  average  length  of  the  two  sides  of  the  column,  and  the 
effective  width  is  the  average  length  of  the  top  and  bottom.  The  resulting  cells  are 
mostly  quadrilaterals,  as  shown  in  Fig.  4.4b.  Some  fill-in  triangular  cells  art,  created 
when  the  number  of  cells  changes  from  one  column  to  another. 


CHAPTER  4.  NUMERICAL  SOLUTIONS  90 


FigTire  4.4:  Domain  Partition,  (a)  First  step  -  column  division;  (b)  second  step  -  cell 
division. 

4.2  Evaluation  of  Integrals 


The  boundary  element  method  depends  on  the  evaluation  of  integral  equations  to 
model  a  system.  Many  types  of  integrands  are  involved,  and  practically  all  of  them 
exhibit  singular  behaviors,  creating  difficulties  for  numerical  implementations.  As  the 
heart  of  the  BEM,  integration  techniques  warrant  a  detailed  examination. 

The  division  of  the  boundary  into  segments  facilitates  the  estimation  of  source 
distribution  and  the  prescription  of  boundary  conditions.  Similarly,  the  representation 
of  segments  by  straight  lines  simplifies  the  evaluation  of  line  integrals.  Even  so,  these 
expressions  are  still  difficult  to  derive,  as  the  line  path  defining  the  integration  limit 
are  arbitrarily  placed  in  a  two-dimensional  space.  The  predicament  becomes  worse  for 
area  integrals  on  quadrilaterals  or  triangular  elements^. 


^Fortunately,  the  pseudo-sources  that  are  used  for  modeling  nonlinear  bulk  parameters  can  be  com- 
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In  two-dimensional  Laplace  problems  that  employ  4>*  (i)  =  log(z)  as  the  kernel  fimc- 
tions,  line  integral  involving  (f>*  and  its  first  derivatives  are  routinely  done  in  analytical 
form.  This  approach  is  usually  abandoned  in  biharmonic  systems  because  the  functions 
are  more  complicated.  For  periodic  problems  that  use  =  log{2(cosh(y)  —  cos(i))] 
kernel,  such  an  anzdytical  evaluation  actually  becomes  infeasible.  The  Lobachevski 
functions  family  L{x)  =  /log[stn(x)]dz,  to  which  the  periodic  kernel  yields  after  in¬ 
tegration,  only  exists  in  integral  form  and  cannot  be  reduced  further  -  just  like  the 
error  function  erf(x).  But  unlike  the  error  function,  no  library  routines  are  available 
to  provide  numerical  results  for  Lobachevski  functions.  Thus  niimerical  integration  is 
used  whenever  possible  and  convenient. 


Consider  the  integral  equations  stated  previously  in  Eqs.  3.7  and  3.11: 


^(p)  =  _^p(q)9^*(p-q)dr 


^^(P)  ^ 
dn 


a<^*(p-q) 

dn 


dr 


After  boundary  segmentation,  they  become 

N 


^(p)  «  /  ^*(p-q)d? 

f=i  ■'r.- 


dn 


1=1 


dn 


dq 


(4.1) 

(4.2) 


where  N  is  the  number  of  segments  and  a,-  is  the  source  distribi  Ion  on  T,-,  the  x** 
boimdajy  segment.  The  a’s  are  factored  out  from  the  integral  expressions  because 
they  are  constants.  Methods  of  evaluating  straight-line  segment  integrals  are  discussed 
in  the  following  sections. 


4.2.1  Numerical  Quadratures 

Numerical  integration,  or  quadrature  as  it  is  often  called,  seeks  to  approximate  the 
integration  //(x)dx  by  evaluating  the  expression  /(x)  at  discrete  points  x,-,  weighing 
puted  by  integrating  around  the  perimeter  of  a  cell,  as  shown  later.  No  area  integrals  are  needed. 
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them  with  Wi  and  summing  up  the  results: 

i-O 

The  two  most  common  schemes  to  determined  the  spacings  of  a:,-  and  values  for  Wi  are 
the  Gaussian  quadrature  and  the  Newton-Cotes  formulae.  In  the  following  discussions, 
it  is  assumed  that  the  integration  goes  from  —1  to  1. 

Newton-Cotes  Formulae 

In  the  Newton-Cotes  formulae,  the  nodes  i,’s  are  equally  spaced.  If  the  endpoints  of 
the  line  path  are  also  quadrature  points,  the  formula  is  said  to  be  closed,  otherwise 
it  is  open.  The  accuracy  or  the  “order”  of  the  formula  depends  on  the  number  of 
quadrature  points.  Of  the  low-order  closed  approximations,  the  2"^  order  is  more 
commonly  known  as  the  trapezoidal  rule,  and  the  3'’'*  order  Simpson’s  rule.  The 
closed  Newton-Cotes  formulae  that  have  been  used  in  this  investigation  are  given  in 
Table  4.1.  Observe  that  the  7*^  order  approximation  has  an  undesirable  property  in 
that  Wi  fluctuates  wildly.  In  fact  the  9*'*  (not  shown)  and  higher  order  approximations 
have  imdesirable  negative  weighing  factors. 

Gaussian  Quadrature 

The  Newton-Cotes  approximation  was  later  replaced  with  Gaussian  quadrature.  The 
Gaussian  sheme  uses  irregulzirly  spaced  nodes.  The  coefficients  for  Gaussian  quadra¬ 
ture  are  given  in  Table  4.2. 

The  reason  for  the  need  of  high  order  approximations  is  that  flelds  have  fast  varying 
fields  near  a  singularity.  This  situation  arises  when  the  top  and  bottom  boundaries 
lie  very  close  together,  as  oxide  is  thin  at  the  beginning  of  oxidation.  The  selection 
of  quadrature  order  is  determined  by  the  closeness  the  collocation  point  is  to  the 
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n 

±Xi 

Wi 

3 

0 

4/3 

1 

1/3 

5 

0 

12/45 

1/2 

32/45 

1 

7/45 

7 

0 

272/420 

1/3 

27/420 

2/3 

216/420 

1 

41/420 

Table  4.1:  Newton-Cotes  Coefficients 


n 

Wi 

■3“ 

0.0000000000 

0.7745966692 

0.8888888888 

0.5555555555 

5 

0.0000000000 

0.5384693101 

0.9061798459 

0.5688888888 

0.4786286704 

0.2369268850 

7 

0.0000000000 

0.4058451513 

0.7415311855 

0.9491079123 

0.4179591836 

0.3818300505 

0.2797053914 

0.1294849661 

9 

0.0000000000 

0.3242534234 

0.6133714327 

0.8360311073 

0.9681602395 

0.3302393550 

0.3123470770 

0.2606106964 

0.1806481606 

0.0812743883 

Table  4.2:  Gaussian  Quadrature  Coefficients 
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boundary  segment®. 

4.2.2  Analytical  Integrations 

Source  contributions  from  all  segments  are  summed  up  to  determine  the  boimdary 
condition  of  a  particular  segment.  When  the  source  coincides  with  tne  collocation 
point,  the  integrand  becomes  indefinite.  Under  such  a  condition,  analytical  integration 
must  be  used  as  numerical  quadratxire  performs  rather  poorly.  If  a  quadrature  node 
coincides  with  the  collocation  point,  the  numerical  solutions  will  blow  up  as  we  attempt 
to  evaluate  log(O)  or  divide  something  by  0.  On  the  other  hand,  if  the  singular  point 
is  omitted,  the  approximation  will  be  in  great  error  because  the  contribution  of  many 
constant  come  from  the  singular  point,  as  we  will  see  later.  Only  analytical  integration 
can  deal  with  singularity  adequately. 

Periodic  kernels,  as  noted  that  in  Chapter  3,  can  be  reduced  to  non-periodic  coun¬ 
terparts  around  the  singularities  by  xising  Taylor  series  expansions  of  the  trigonomet¬ 
ric  and  hyperbolic  functions.  All  analytical  integrations  over  singularities  are  done 
in  terms  of  those  non-periodic  versions.  Because  of  the  limitations  of  Taylor  series 
approximation,  the  integral  must  be  done  close  to  the  singular  point,  and  the  path 
cannot  be  too  long.  Typically  the  segment  length  is  limited  to  be  shorter  than  0.2 
radian.  For  a  longer  boundary  segment.  It  is  possible  to  use  quadrature  on  the  edges 
and  analytical  integration  at  the  center. 


Coordinate  Transformation 


In  numerical  quadrature,  it  does  not  matter  in  what  order  the  contributions  from 

the  quadrature  points  are  summed,  the  results  are  the  same.  Such  is  not  the  case 

^Not«  that  only  odd  nainb«r  of  nodes  are  employed  m  either  Newton-Cotes  or  Gaussian  quadrature. 
This  has  to  do  with  the  programming  feature  that  checks  the  central  node  before  deciding  which  order 
of  quadrature  should  be  used. 
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for  analytical  integration  -  the  solution  depends  on  the  direction  of  integration  path. 
Incorporating  all  these  information  will  make  the  resulting  expressions  unnecessar¬ 
ily  complicated.  Thxis  a  modular  approach  is  chosen  to  breakdown  the  integration 
procedure  into  two  parts.  First,  the  integral  expressions  are  solved  in  a  local  x'  —  y' 
coordinate  system.  This  local  frame  is  chosen  such  that  the  integration  path  always 
runs  in  a  predefined  direction  of  either  s'  or  j/.  The  solutions  are  then  transformed 


back  to  the  global  reference  frame  using  the  following  rules: 

X  =  cos  dx'  —  sin  Oy*  (4.3) 

y  =  sin  ^  x' +  cos  (4.4) 

X*  =  cos*  dx^  —  cos  B  sin  Bx'y'  +  sin*  By'^  (4.5) 

xy  =  cos*  tfx'y'  +  cos  B  sin  fl(x'*  —  y'*)  —  sin*  5x'y'  (4.6) 

y*  =  cos*  dy'*  +  cos  B  sin  Bx'y'  +  sin*  l9x'*  (4.7) 

X*  =  cos*  5 x'*  —  3  cos*  ^  sin  5x'*y'  4-  cos  B  sin*  x'y**  —  sin*  fly**  (4.8) 

x*y  =  cos*  Bx'^y'  +  cos*  B  sin  (?(x'*  —  2x'y^)  +  cos  8  sin*  6(y'^  —  2x'*y') 

-t- sin*  ^x'y'*  (4.9) 

xy*  =  cos*  dx'y'*  +  cos*  5  sin5(2x'*y' —  y**)  +  cosffsin*  fl(x'*  —  2x'y'*) 

—  sin*  5x'*y'  (4-10) 

y*  =  cos*  0y'*  +  3cos*  ^  sintfx'y'*  +  3cos0sin*  5x'*y' +  sin*  Ax'*  {4-11) 


where  B  is  the  rotation  angle  from  the  x  axis  to  the  x*  axis,  measured  counterclock¬ 
wise.  Note  that  these  rules  also  applies  to  differentiation  operators  by  substituting 
parameters  with  differentiation  operators.  For  instance: 


— T  =  cos'*  B~^  ~  cos  B  sm  B 
ax^  ar* 


.  2 

5x^53/ ^dy'^' 


Working  in  a  local  frame  where  the  integration  path  runs  either  in  the  x'  or  y'  direction 
greatly  simplifies  the  anzdytical  expressions. 
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4.2.3  Line  Integrals  Over  Collocation  Points 

It  turns  out  that  the  integrals  over  the  collocation  point  usually  yield  simple  results 
due  to  the  midpoint  collocation  method  and  the  step  fxmction  approximation  of  the 
sources.  Great  care  must  be  taken,  however,  when  dealing  with  singularities.  The 
integral  path  must  be  chosen  carefully  to  be  consistent  with  the  source-free  assumption 
in  the  boundzuy  integral  equations. 


Three  types  of  integrands  are  encountered.  First,  for  non-singular  or  weak  sing\ilar 
functions  such  as  logarithm  log(r),  we  have  a  well-defined  and  well-behaved  solution. 
For  strong  singularity,  we  get  two  other  types  of  behaviors.  The  second  one  comes 
from  expressions  of  the  form 

- rdx  =  /  -di 

■a  x’  +  y*  J-A  X 

=  [log{x)]t^ 

=  0. 


Although  the  overall  result  is  zero,  it  Is  due  to  cancellation  of  two  logarithm  terms 
that  become  singular  as  a  — *■  0.  This  integration  must  be  interpreted  in  the  Cauchy 
principle  value  sense.  As  we  can  see  from  this  expression,  the  result  shows  logarithmic 
singular  behavior  if  one  of  the  endpoints  is  0  instead  of  ±a.  This  is  the  cause  of  singular 
behaviors  when  the  surface  orientation  is  discontinuous.  The  third  integration  type  is 
of  the  form 


/: 


-dx;  y  -*■  0, 


I -a  x’  -h  y* 

which  yields  0,  tt,  or  —tt,  depending  y  =  0,0'^,  or  0“.  We  must  choose  y  to  be  O'*"  or  0~ 


such  that  the  source  at  the  singularity  jiist  lies  ‘outside’  the  simulation  region  B.  We 
don’t  want  any  source  inside  the  region’.  For  this  expression,  the  jump  in  values  is 

^Actually,  y  is  chosen  to  be  exactly  sero  in  many  boundary  element  references  [38,42,40];  the  integral 
is  0.  But  now  the  half  of  the  source  point  is  lying  inside  the  domain,  therefore  the  r  term  is  picked  up 
by  the  divergence  operate  on  the  source.  We  see  that  these  two  approaches  are  consistent  with  one 
another. 
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due  to  the  term  tan"^  (i/y),  which  is  problematic  in  many  applications,  as  it  suddenly 
changes  values. 

The  integral  expressions  for  the  basis  functions  defined  in  Chapter  3  are  given  in 
the  Appendix. 

4.2.4  General  Analytical  Line  Integrals 

Over  the  course  of  investigating  different  methods  of  modeling  nonlinearity,  it  was 
found  that  high  accuracy  in  the  integral  technique  is  desirable,  not  to  avoid  numerical 
errors,  but  to  enstire  that  integration  errors  do  not  play  a  role  in  the  ill  behavior  of 
the  solutions.  A  set  of  analytical  integration  formulae  is  also  given  in  the  Appendix. 
However  they  are  only  applicable  to  the  non-periodic  basis  functions. 


4.3  Matrix  solutions 

In  Eqs.  4. 1-4.2,  a’s  are  the  unknowns.  Values  for  cr’s  must  be  calculated  before  the 
domain  solution  or  the  unknown  boundary  parameters  can  be  determined.  The  pro¬ 
cedure  is  to  compute  o’s  in  terms  of  the  “known”  boundary  conditions.  Suppose  we 
axe  dealing  with  a  Dirichlet  problem  in  which  $  at  the  boundary  is  specified.  the 
potential  at  the  collocation  point  p,-  of  segment  t,  is  given  as 

)  =  - 
j=i 

= 

where  a,-,-  ;=  —  ?)dr.  Combining  this  expression  from  all  boundary  elements, 

the  following  matrix  representation  is  arrived: 


Ax  =  6. 


(4.12) 
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Ais  3.  N  X  N  matrix  whose  entry  Oj,-  is  defined  above,  and  x  and  b  are  vectors  whose 
elements  axe  {<ri  •  •  •  Cff}  and  •  •  •  $jsr  respectively.  The  solution  of  x  is  given  by 

X  =  (4.13) 

With  cr’s  solved,  we  can  proceed  to  calculate  at  the  boundary. 

Suppose  the  problem  is  Neumann,  Eq.  4.2  is  used  to  assemble  a  similar  matrix 
equation  to  solve  for  a’s: 

Cx  =  d.  (4.14) 

where  c,-,-  =  5<^*(pi  -  q) /-^ndq,  and  d  =  d^i/dn,  •  •  • ,  5$Ar/5n.  For  simplicity,  we 

have  ignored  the  required  normalization  equation. 

4.3.1  Matrix  System  for  Diffusion 

The  oxidant  diffusion  problem  is  a  mixed  boundary-value  problem  consisting  of  Dirich- 
let,  Neumann  and  Robin  botmdary  condition,  as  mentioned  in  Chapter  2.  A  new 
hybrid  matrix  equation  is  assembled  from  Eqs.  4.12  and  4.14,  given  symbolically  as 

(A|C1x  =  615.  (4.15) 

Depending  on  the  boundary  condition  of  a  particular  segment  t,  the  t**  row  of  A  or  C, 
and  element  of  b  or  d,  or  a  linear  combination  of  both,  is  used  to  form  the  entries 
of  this  matrix: 

ttii  • • •  Cis 

•  •  • 

«  %  • 

•  •  • 

Cii  •  •  •  Ciif 

•  •  • 

•  •  • 

•  •  • 

where,  for  illustration,  boundary  elements  1  and  I  are  shown  to  have  a  Dirichlet  con¬ 
dition  and  a  Neumann  condition  respectively. 

Note  that  all  those  matrices  A,  C  and  the  hybrid,  are  fully  populated.  This  is  due 
to  the  global  influence  of  the  kernels.  The  somce  a  on  any  one  segment  has  a  direct 
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effect  on  the  boundary  condition  of  on  all  segments.  In  contrast,  the  finite  difference 
and  finite  element  methods  always  produce  sparsely  populated  matrix  because  of  the 
short-ranged  interaction  between  elements.  Their  nodes  or  elements  only  interact  with 
nearest  neighbors. 

The  solution  for  the  imknown  x  is  obtained  using  Gaussian  or  Gauss-Jordan  elimi¬ 
nation  with  partial  pivoting.  Note  that  an  explicit  inversion  of  the  matrices  is  seldom 
required, 

4.3.2  Oxide  Flow 

The  boundary  conditions  for  oxide  flow  are  specified  in  vector  form.  For  a  two- 
dimensional  problem,  each  boundary  segment  requires  two  entries  in  the  matrix.  Thus, 
the  resulting  matrix  size  is  2N  x  2N.  Otherwise,  the  matrix  synthesis  process  is  similar 
to  that  of  diffusion. 

4.3.3  Multiple  Domains 

The  discussion  so  far  concerns  with  single-domain  problems.  This  treatment  is  ade¬ 
quate  for  oxidant  diffusion  but  not  for  the  oxide  motion.  The  oxide  bulk  is  in  contact 
with  silicon  nitride  mask  and  silicon  substrate  (and  also  the  ambient).  For  the  diffusion 
problem,  the  boundary  conditions  at  these  interfaces  are  F  •  ri  =  0  autid  kgC  =  DF  •  n, 
due  to  the  fact  that  silicon  nitride  is  impermeable  to  oxidant  and  that  chemical  re¬ 
action  is  a  interfacial  activity.  The  mechamical  interaction  between  these  bodies  are 
more  involved.  All  materials  are  flexible.  But  to  medse  the  problem  simpler,  we  as¬ 
sume  that  the  silicon  substrate,  which  is  very  thick  with  respect  to  the  oxide,  is  totally 
rigid.  However,  the  silicon  nitride  mask  is  of  comparable  thickness  to  oxide  and  de¬ 
forms  dxiring  oxidation,  therefore  its  interaction  with  the  oxide  motion  may  not  be 
overlooked. 
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At  the  interface  between  oxide  and  nitride,  neither  the  displacement  nor  surface 
traction  can  be  preassigned.  The  actual  values  for  these  parameters  depends  on  the  de¬ 
tail  interactions  of  these  bodies.  A  general  approach  for  dealing  with  multiple  domains 
is  to  combine  their  matrix  equations  into  one  amd  match  the  boundary  conditions  at 
the  common  interface.  Let’s  assume  that  the  matrix  equations  for  two  domains,  A 
and  B,  are  given  as 

\Gf' 

[Hi- 

[iff. 

where  C/’s,  P’s,  and  o’s  are  the  surface  displacement,  traction,  and  sources  respectively. 
The  domains  are  identified  by  the  superscripts;  the  boundaries  are  denoted  by  the 
subscripts. 


(-) 

^  =  (  yC  I  (4.19) 

(4.20) 


The  requirements  at  the  common  interface  T/  are  (i)  displacements  be  continuous: 

Uf  =  Uf 

and  (ii)  surface  tractions  be  equal  and  opposite: 


The  first  ensures  no  overlapping  of  materials  is  possible,  the  second  conserves  forces. 
Note  that  these  vectors  are  defined  with  respect  to  the  global  reference  frame  for 
consistency,  instead  of  the  local  frames.  For  the  combined  systems,  we  thus  have  the 
following  matrbc: 


[G\H]t  0 

\G]f  -[G]l 

mf  [ml 

0  (G1  J]f 


0 

0 

V  {[u\P]}^ 


(4.21) 
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This  matrix  is  larger  but  contains  many  zero  entries.  It  has  been  found  that  the  matrix 
solution  can  be  sped  up  significantly  by  taking  advantage  of  the  blocked  nature  of  the 
matrix. 

4.3.4  General  Domain  contributions 

In  the  interior-cell  formulation  for  nonlinear  diffusion  and  oxide  motion,  there  is  no 
difference  between  the  contributions  from  the  domain  sources  and  those  from  the 
bovmdary  somces,  as  far  as  the  matrix  setup  is  concerned.  For  an  equivalent  system, 
the  interior-cell  approach  produces  a  denser  and  less  efficient  matrix  thzm  the  subregion 
technique,  however  it  is  more  readily  extended. 


4,4  Computation  for  Nonhomogeneous  Domain 


By  inducing  variations  in  the  diffiisivity  and  viscosity,  stress  causes  pseudo  domain- 
sources  to  appear  in  the  integral  equations  for  both  oxidant  diffusion  and  oxide  motion 
processes.  Simplifications  and  techniques  for  evaluating  the  domain  contribution  axe 
outlined  below. 

We  have  mentioned  in  Chapter  3  that  we  use  daf^/dxj  and  V  •  F  instead  of  the 
nonderivative  forms  in  our  domain  calculations.  The  reason  has  to  do  with  the  domain 
integral  approximation.  In  our  scheme,  the  domain  is  partitioned  into  triangular  or 
quadrilateral  cells,  as  described  earlier  in  this  chapter.  Within  each  cell,  it  is  assumed 
that  the  scale  factors  are  constant;  they  take  finite  jumps  from  one  cell  to  another. 
Consequently  the  values  for  those  pseudo  sources  are  zero  inside  a  cell  but  become  an 
impulse  at  the  cell  boundaries.  This  features  reduces  domain  area  integrals  into  line 
integrals  along  the  perimeter  of  the  cells. 

The  values  for  the  scale  factors  are  computed  using  stress  values  at  the  center  of 
the  cell.  For  the  source  strength  of  the  cells,  two  different  estimation  techniques  have 
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© 
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(b) 


Figure  4.5:  Dipole  sources.  Unlike  a  simple  source,  a  dipole  source  can  have  spatial 
orientation. 

been  considered. 

4.4.1  Dipole  and  Couplet  Approximations 

As  a  vehicle  to  realize  nonlinearities  in  the  bulk  parameters,  the  pseudo  domain-sources 
are  supposed  to  perturb  the  solution,  but  not  to  introduce  real  body-force  or  oxidant 
source  in  the  bulk.  Hence  the  net  values  of  the  domain  sources  should  be  zero. 

One  way  to  automatically  satisfy  the  zero  average-value  constraint  is  to  employ 
dipole  and  couplets.  As  the  kernel  fimction  for  the  double-layer  method  mentioned  in 
Chapter  3,  a  dipole  is  obtained  from  <^*  by  taking  derivative  in  an  appropriate  direction. 
Shown  in  Fig.  fg:dipoles  are  the  two  basic  dipoles,  aligned  in  the  x  and  y  directions, 
respectively.  Because  positive  and  negative  charges  (source  and  sink)  always  appear 
in  pair,  the  flux  contribution  to  an  enclosing  boundary  is  zero.  The  couplets  are  a 
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(a)  (b)  (c) 


Figure  4.6:  Couplet  sources 

pair  of  force  acting  in  opposite  directions.  The  three  basic  couplet  configurations  are 
shown  in  Fig.  4.6. 

When  one  dipole  or  couplet  soiirce  is  distributed  uniformly  over  a  domain  fl,  the 
positive  and  negative  charges  overlaps  and  cancels  in  the  bulk,  leaving  unbalanced  ones 
on  the  perimeter.  This  situation  is  shown  graphically  in  Fig.  4.7  for  a  dipole  source. 
It  can  be  shown  mathematically  that,  with  integration  by  pzirts,  an  area  integral  is 
reduced  to  a  line  integral  around  the  perimeter. 

In  the  first  approach  for  determining  pseudo  domain-source,  the  flux  and  stress 
are  estimated  at  the  cell  center  and  assumed  to  be  uniform  across  the  cell,  just  like 
the  scale  factors  and  Each  cell  is  considered  individually  without  regards  to 
neighbor  cells.  The  resulting  source  distributions  resemble  the  dipoles  and  couplets, 
thus  the  name  for  for  this  correct  method. 


When  used  on  simple  test  structures  where  diffusion  flux  variation  is  smooth  {  for 
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Domain  Source  Boundary  S  urface  Source 


Figure  4.7:  Conversion  of  Dipole  Area  Integral  into  Contotir  Integral.  The  transfor¬ 
mation  is  possible  due  to  the  fact  the  dipoles  in  the  interior  cancel  one  another. 

the  homogeneous  case),  the  performance  of  this  method  is  adequate,  t.e.  it  can  block 
the  diffusion  to  a  good  extent  if  the  scale  factors  kd's  for  the  test  cells  are  reduced 
significantly.  Some  residual  fiux  is  present  unavoidably  as  this  is  a  perturbation  tech¬ 
nique.  Unfortunately  for  real  oxidation  structures  such  as  LOCOS,  the  performance  is 
rather  dismal.  What  happens  is  that  the  oxidant  fiux  near  the  edge  of  the  nitride  mask 
is  highly  nonuniform.  Because  of  the  thinness  of  the  relief  oxide,  interior  cells  there 
are  only  single-layer.  Sources  determined  at  the  cell  center  cannot  provided  adequate 
correction  to  the  flux  distribution  at  the  cell  boundaries  -  it  is  not  sufficient  to  coun¬ 
teract  the  oxidant  flux  reaching  the  silicon  side,  jmd  on  the  other  hand,  oversupply 
oxidant  to  the  nitride  side,  which  is  reflected  back.  Consequently,  it  losses  its  ability 
to  block  the  diffusion.  In  fact,  it  produces  terrible  results  -  where  the  cell  partition 
ends,  the  supply  of  oxidant  may  increase.  The  performance  of  the  doublets  on  the 
oxide  flow  problem  is  equally  lackluster  --  it  czm  never  reduce  the  stress  level  by  more 
than  half  regeu'dless  how  small  the  scale  faotor  is. 
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4.4.2  Interfacial  Source  Formulation 

One  major  weakness  with  the  dipole/couplet  scheme  is  that  the  collocation  point  for 
the  domain  sources  is  located  at  the  cell  center  -  flux  and  stress  are  assumed  to 
be  uniform  across  the  cell.  This  problem  is  circumvented  in  the  interfacial  source 
formulation.  Major  activities  such  as  oxidation  take  places  on  the  sides  of  a  cell.  By 
having  the  source  computed  individually  at  each  side  of  a  cell,  wider  variations  of 
source  strength  can  be  achieved  to  cope  with  the  highly  nonuniform  flux. 

To  decide  how  the  interfacial  source  is  computed,  we  consider  the  domain  integrzd 
of  Eq.  3.93: 

At  the  interface  between  two  cells,  the  value  of  kd  changes  abruptly  from  ki  to  k2, 
resulting  in  non-conservation  of  flux.  Given  the  fact  the  only  the  normal  component 
of  flux  matters  and  that  the  interface  is  a  straight  line,  at  locations  very  close  to  the 
interface,  the  behavior  is  one  dimension?!. 

Hence  we  examine  an  equivalent  one-dimensional  model,  -using  physical  arguments 
instead  of  mathematics.  As  shown  in  Fig.  4.8,  this  model  has  two  regions  -  1  and  2  - 
having  different  diffusivity  ki  and  ki*  (For  simplicity,  we  drop  Dq  as  it  never  appears 
in  the  flnal  solution.)  The  flux  In  each  region  is  given  by 

Fi  =  kiC 


and 


Fz  =  kzC 


respectively.  Thus,  at  the  interface,  there  is  a  build-up  of  oxidant  of 


{ki-k2)C. 


which  has  to  be  dissipated.  Green’s  function  for  one-dimensioned  diffusion  is  —  |z|, 
which  has  a  gradient  of  ±1,  depending  on  whether  x  is  smaller  or  larger  than  zero. 
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Figure  4.8:  One-Dimensional  Nonhomogeneous  Diffusion.  The  secondary  fluxes  are 
created  due  to  the  mismatch  in  diffusivity  between  the  two  regions. 


Hence  the  fluxes  generated  by  the  interface  are 

FI  =  -k,D 
F\  =  k^D 


where  D  is  the  source  strength  yet  to  be  determined  and  the  minus  sign  in  FI  indicates 
that  the  flux  goes  in  the  negative  direction,  as  illustrated  at  the  bottom  of  Fig.  4.8. 
Therefore  net  outgoing  flux  is  given  by  (Fj  —  FI),  and  it  must  equal  to  net  incoming 
flux,  which  is  (Fi  —  F2).  Thus  we  have  the  following  expression  for  the  interfacial 
source: 


D^hzhc. 


(4.22) 


ki  +  *2 

Simply  speaking,  when  a  flxix  impedes  on  a  barrier  (that  has  lower  diffusivity),  it  goes 
through  partially.  Due  to  build-up  of  materials  at  the  boundairy,  two  additional  fluxes 
are  created:  one  is  reflected  back  (F^),  the  other  continues  forwards  (Fj).  These  fluxes 
seem  to  be  created  by  a  pseudo  domain-source.  Note  that  even  if  fc:  is  zero,  t.e.  when 
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region  2  is  impermeable,  the  oxidant  concentration  in  it  is  nonzero  and  show  a  finite 
gradient.  This  is  due  to  the  steady-state  difFtision  assumption.  This  technique  has 
been  verified  to  work  in  one  dimension  analytically,  and  two  dimensions,  numerically. 

The  extension  to  viscoelastic  flow  is  straight  forward.  The  effective  stress  scale 
factor  is  and  kt,2  for  Region  1  and  2  respectively.  The  tinscaled  stress  distribution 
is  continuotis  across  the  interface.  Each  region  contributes  towards  the  surface  traction 
on  the  interface.  If  k„i  =  kg^,  the  contributions  will  cancel  each  other.  On  the  other 
hand,  if  kgi  kg^,  3,  net  force  will  be  exerted  on  the  interface.  It  is  nullified  by  a 
pseudo  body-force  which  generates  its  own  stress  and  displacement  fields.  The  value 
of  the  body-force  vector  is  given  by 

kgl  -  kg2 

T  I 

Kgl  -p  Kg2 

As  we  can  see,  the  nonlinear  viscoelasticity  formulation  parallels  that  of  nonlinezir 
diffusion. 

A  last  note,  the  calculation  for  a  particular  pseudo  source  must  not  include  the 
flux  or  stress  generated  by  the  source  itself.  This  is  can  be  done  easily  in  the  program 
codes. 

Cells  in  Contact  with  Boundary  Segments 

With  the  formulae  given  by  Eqs.  4.22  and  4.23,  we  can  calculate  the  interfacial  source 
between  two  cells  and  between  a  cell  and  a  homogeneous  region.  One  remaining  ques¬ 
tion  is  on  the  proper  treatment  for  a  cell  connected  to  a  boundary  segment.  Recall 
the  boimdary  source  is  located  on  the  outside  of  a  boundary  collocation.  The  inter- 
facial  source  is  by  definition  on  the  inside.  Shown  schematically  in  Fig.  4.9a.  are  the 
positions  of  the  sources  with  respect  to  the  collocation  point.  The  distances  ,  Si  and 
S2,  are  infinitesimally  small.  Consider  the  diffiision  problem.  The  sandwich  area  con¬ 
taining  the  collocation  point  is  outside  a  cell,  and  therefore  assumed  to  belong  to  the 
homogeneous  confine.  Its  diffusivity  is  the  nominal  value  D  =  Do,  but  this  really  does 
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Figure  4.9:  Boundary  Sources. 


(a)  Original;  (b)  merged. 


not  matter  because  the  layer  is  so  thin  that  the  transport  process  is  dominated  by  the 
cell.  Given  the  close  proximity  of  the  two  sources  (di  +  5:  — »  0),  one  would  expect 
that  the  value  of  one  source  would  track  the  other  in  some  specific  way.  Thtis  it  is 
useful  and  beneficial  to  recode  the  mathematics  so  that  one  source  (t.e.  the  interfacial 
one)  can  be  eliminated.  Through  derivations  that  au’e  not  presented  here,  indeed,  the 
interfacial  source  cam  be  eliminated  and  in  place,  and  assign  the  boundary  segment 
the  effective  diffusivity  D  =  Dq.  The  new  configuration  is  shown  in  Fig.  4.9b.  The 
same  approach  works  for  the  nonlinear  viscoelasticity  problem. 


4.5  Nonlinear  Solutions 

The  effects  of  stress  on  thermal  oxidation  may  be  divided  into  two  categories,  depend¬ 
ing  on  whether  they  act  on  the  boundary  conditions  or  on  the  bulk  parameters.  Either 
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case,  they  are  nonlinear  functions  that  have  to  be  solved  iteratively.  The  two  boundary 
parameters  influenced  by  stress  are  the  surface  reaction  rate  kg  and  the  equilibrium 
concentration  C*  of  the  diffusion  system.  Because  boxmdary  conditions  do  not  directly 
involve  the  domain  flux  conservation  law,  the  system  remains  Laplacian  and  is  still 
subject  to  reduced  dimensionality  treatment.  Bulk  parameters  affected  by  stress  in¬ 
clude  the  diffusion  coeSicient  of  oxidants  and  the  viscosity  of  oxide.  For  them,  domain 
corrections  outlined  in  the  previoiis  section  are  used. 

The  stress  functions  used  in  this  investigation  are  of  this  form: 


where  P  is  a  specific  component  of  the  stress  tensor  ,  V  is  an  activation  volume,  Jfc 
is  Boltzmann’s  constant,  and  T  is  the  temperature.  In  descriptions  that  follow,  it  is 
assumed  that  the  system  is  stable  (t.e.  it  does  not  have  positive  feedback)  and  that 
there  is  only  one  solution. 

4.5.1  Scalar  Iterative  Methods 

All  the  stress  effects  requires  the  solution  to  a  nonlinear  function,  conceptually  of  the 
following  form: 

/(x)  =  X  (4.24) 

where,  for  a  simple  case,  x  and  /  are  scalars.  A  few  iterative  techniques  are  available 
for  obtained  the  solution  for  x.  One  basic  approach  is  to  use  this  updating  scheme: 

x‘+'  =  /(x‘). 

where  the  superscripts  indicate  the  iteration  numbers.  x°  is  a  starting  guess  value. 
This  technique  shows  poor  convergence,  nonetheless  it  is  simple  and  does  not  require 
additional  information  of  the  system. 
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Newton’s  Method 


For  faster  and  more  stable  convergence,  Newton’s  method  is  often  employed.  To 
illustrate  this  technique,  we  recast  Eq.  4.24  in  a  dilferent  form: 

g(x)  =  /(x)  -  X. 

where  j(x)  is  the  residue.  g(x)  is  0  if  x  is  the  root  to  Eq.  4.24.  Suppose  Xo  is  the 
solution  (which  we  don’t  know  yet),  and  x  is  close  to  xo,  taking  the  Taylor  series 
expansion  of  ff(xo)  to  the  1**  order  at  x  yields 

gixo)  «  p(x)  +  p'(x)(xo-x)  (4.25) 

«  0. 

where  the  prime  (')  denotes  the  first  derivative  with  respect  to  x.  Accordingly,  xq  as 
X  —  g(x)/g'(x).  We  have  thus  derived  the  Newton  iteration  scheme: 

=:  X*  -  (^-26) 

p’(x*) 

This  scheme  shows  a  quadratic  convergence,  i.e.  error  is  reduced  by  half  in  every 
iteration.  However,  it  still  may  fail  to  converge  under  some  conditions.  For  instance, 
the  solution  may  oscillate  between  two  regions  for  certain  /(x)’s  and  initial  conditions. 
Also  the  solution  can  blow  up  if  g'(x)  happens  to  be  zero. 


It  is  difficult  to  come  up  with  a  general  iterative  method  that  is  always  robust 
and  converges  quickly  for  any  type  of  /(x).  Although  not  always  possible,  one  should 
analyze  the  problem  nature  and  determine  the  best  possible  solution  approach.  Indeed 
the  main  reason  that  the  Newton  scheme  has  a  fast  convergence  is  that  it  considers 
the  direction  the  system  is  heading,  namely  /'(x).  In  case  /'(x)  cannot  be  obtained 
directly,  one  may  approximate  it  by 


/'(x‘) 


_  /(x*)-/(x^ 


-1^ 


(4.27) 


The  iteration  scheme  thus  become: 


X*  —  x*~^ 
5(x*)  -  p(x*-i) 


i 


x*"^^  =  X*  —  ^(x*) 


(4.28) 
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which  is  known  as  the  secant  or  false  position  technique  [44]. 


Modified  Iterations 


Sometimes  it  is  necessary  to  "damp”  the  iterations  in  order  to  stabilize  the  solution. 
This  is  achieved  by  adding  a  relaxation  parameter  (O  <  A  <  l)  to  scale  down  the 
change  in  x  at  each  iteration.  The  modification  is 


=  (1  -  A)2‘  +  A/(x‘) 
for  the  simple  updating  scheme  and 

=  X*  -  Ag(x*) 

j(xfc)  -  g(x*-i) 

for  Newton’s  method. 


(4.29) 


(4.30) 


4.5.2  Multi- Variable  Solutions 

The  scalar  iteration  methods  can  be  extended  to  desd  with  multiple  variables,  but 
there  are  complications.  For  one  thing,  the  variables  in  general  are  not  decoupled.  In 
other  words,  changing  the  condition  on  one  variable  affects  the  rest.  It  is  difficult  to 
predict  in  what  direction  the  the  system  is  moving.  Consequently  schemes  such  as  the 
simple  updating  may  easily  become  unstable  if  they  update  the  solution  in  the  wrong 
direction.  The  solution  may  need  to  be  damped,  resulting  in  slow  convergence. 


Updating  Scheme 

With  the  damping  factor  A  incorporated,  the  multi-variable  updating  scheme  is 

x*-*-^  =  (1  -  A‘)x*  +  A‘7(x*).  (4.31) 

which  is  a  straightforward  extension  of  the  scalar  case. 
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Newton’s  Method 

For  Newton’s  method,  the  vector  Taylor  series  expansion  becomes 

g(xo)  «  Sf(f)  +  J(z)(zi)  -  i)  (4.32) 

«  0. 

where  g(x)  =  /(x)  —  x,  and  J  is  the  Jacobian  of  g. 

Assiiming  that  the  vector  size  of  g  and  xis  N,  the  definition  for  J  is 

dgi 

-STi  ^ 

:  i  (4.33) 

ofiT 

where  the  element  number  is  given  by  the  subscript.  The  Jacobian  is  a  measure  of 
how  g  will  respond  to  a  small  change  in  x. 

Newton’s  method  thus  becomes 

=  X*  -  (i*)^ (z*)  (4.34) 

The  detenmnation  of  the  Jacobian  is  not  trivial  -  the  array  is  large  [N  x  N)  and 
its  components  are  not  readily  obtained  from  g  in  analytical  form.  Using  the  false 
position  technique  is  almost  out  of  the  question  because  enormous  computation  is 
required;  only  rough  approximation  is  used  here.  Even  when  an  accurate  Jacobian 
is  available,  the  method  often  fails  converge.  One  often-cited  reason  is  that  it  tends 
to  overshoot  in  correction  [45].  The  Taylor  series  expansion  is  good  only  when  x  is 
very  close  to  the  actual  solution.  If  it  is  not  close,  the  Jacobian  produces  a  poor 
correction  term.  In  our  applications,  an  error  can  result  in  larger  deviations  due  to 
the  exponential  terms  embedded  in  g. 

Newton’s  scheme  is  modified,  as  suggested  in  [45] : 

+ jj/j-'jtf*) 


(4.35) 
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to  mclude  two  parameters  —  3*  and  tj,  —  to  increase  damping.  The  function  of  tk  is 
the  same  as  A*  in  the  scalar  case  -  to  reduce  the  changes.  Sk  makes  the  Jacobiaui  more 
diagonally  dominant,  effectively  reducing  the  interaction  between  pau’ameters. 

In  both  iteration  techniques,  there  is  no  algorithm  for  calculating  the  optimum 
values  for  the  relaxation  parameters.  They  axe  determined  on  a  triad  and  error  basis. 
A  two-level  iteration  scheme  is  adopted.  The  purpose  of  having  the  sub-iteration  is  to 
find  large  relaxation  parameters  that  does  not  cause  the  system  to  go  unstable. 


4.5.3  Determination  of  Relaxation  Parameters 


There  is  no  formula  to  determine  optimum  values  for  the  relaxation  parameters  A*,  s*, 
and  a  priori.  They  are  obtained  on  a  trial-and-error  basis  at  each  and  iteration.  If 
these  parameters  are  too  small,  then  the  speed  of  convergence  will  suffer.  Conversely, 
if  they  are  too  large,  the  solutions  will  become  unstable. 

These  parameters  are  adjusted  according  to  the  convergence  rate  of  the  solutions. 
Two  types  of  error  defiiution  are  employed  in  our  iteration  for  control  purposes;  they 
are  the  maximum  norm  (||  •  ||oa)  and  the  normalized  Euclidean  norm  (||  •  Us)  defined 
below: 


llflilco 

mu 


maxlyil 


Vn\ 


Zal 

»=i 


(4.36) 

(4.37) 


Note  that  our  definition  for  the  Euclidean  norm  has  been  normalized  by  IfVN  so 
that  it  yields  the  same  result  as  the  maximum  norm  if  all  elements  of  g  are  the  same 
magnitude. 


The  maximum  norm  is  applied  to  a  residual  error  vector  g  to  to  determine  whether 
an  iteration  loop  should  be  terminated  because  the  solution  has  converged,  t.e.  the 
maximum  error  has  dropped  below  a  certain  threshold. 
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The  Euclidean  norm  is  used  for  determining  the  relaxation  parameters.  Since 
it  weights  all  elements  in  g,  the  Euclidean  norm  is  less  sensitive  to  fluctuations  of 
individual  elements  than  the  maximum  norm.  This  feature  is  higher  desirable  because 
in  the  iteration  process,  some  errors  may  arise  locally  and  momentarily. 

The  value  for  a  relaxation  parameter  is  adjixsted  by  comparing  and  ||^|[2 

where  *  denotes  the  current  trial  step  for  iteration  k.  If  the  current  trial  error  is  smaller, 
X*  is  made  and  the  relaxation  parameter  is  increased  slightly  for  the  next  iteration. 
On  the  other  hand,  if  the  current  trial  error  is  larger,  the  relaxation  parameter  is 
reduced  to  increase  damping,  and  x*  and  everything  else  is  recomputed.  This  sub¬ 
iteration  is  repeated  until  the  ctirrent  trial  error  becomes  smaller.  The  sub-iteration 
is  also  terminated  after  a  few  trials  if  it  seems  impossible  to  make  the  error  smaller  at 
the  current  iteration. 

4.5.4  Iteration  on  Reaction  Rate 

The  stress  effect  on  reaction  rate,  A;«,  is  a  self-regulating,  negative  feedback  system  that 
tends  to  make  oxide  growth  rate  more  uniform  along  the  silicon  interface.  Nonlinear  k, 
is  a  surprisingly  difficult  problem  to  solve  numerically.  With  simple  iterative  schemes, 
a  small  local  error  can  easily  transform  into  a  global  instability. 

Let’s  suppose  the  value  of  k,  of  a  particular  boundary  segment  is  higher  than 
its  equilibrium.  Being  more  reactive,  the  segment  therefore  generates  more  oxide 
than  its  two  neighbors;  it  also  deplete  to  a  certain  degree  the  supply  of  oxidants  to 
them.  (We  ignore  segments  further  away.)  Because  it  injects  more  oxide  into  the 
b\ilk,  the  segment  therefore  experiences  a  compressive  surface  traction  that  reduces 
its  kf.  Meanwhile  the  slower-moving  oxide  of  the  two  neighboring  segments  is  being 
pulled,  procducing  a  tensile  surface  traction  on  the  segment.  In  the  next  iteration 
step  of  a  simple  updating  scheme,  the  k,  of  the  segment  in  question  will  be  lowered, 
while  those  of  the  two  neighbors  will  be  raised:  a  reversal  process  is  in  progress.  As 
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iteration  proceeds,  a  ripple  of  alternating  stress  and  will  spread  through  all  the 
segments.  If  damping  is  insufiBcient,  the  errors  become  a  self-reinforcing  wave  that 
grows  in  magnitude,  resiilting  in  global  oscillation.  For  large  damping  factors,  the 
errors  will  still  propagate  to  adjacent  segments,  but  with  decaying  magnitude;  the 
solution  eventually  converges  after  many  iterations. 

The  reason  for  the  poor  performance  is  that  the  technique  does  not  consider  the 
behavior  of  the  system.  In  particular  it  does  not  anticipate  that  a  change  in  ifc,  of 
a  segment  will  affect  others  as  well.  Better  stability  is  only  ensured  when  the  sys¬ 
tem  response  is  factored  into  the  updating  scheme.  This  brings  us  back  to  Newton’s 
method. 

Accordingly,  /  for  normalized  kg  is  of  the  form: 

/  = 

where  p,-  •  A,-  is  the  normal  surface  traction  at  boundary  segment  t.  As  shown,  f  has 
entries  for  M  silicon  interface  segments.  The  surface  traction  list  p  =  {pi,P2,  •  •  •  ,PAf} 
may  be  explicitly  defined  as  a  function  of  the  normalized  Jfc,  list:  p  =  p(7).  As 
a  combined  product  of  the  oxidant  diffusion  and  oxide  fiow  systems,  the  feedback 
mechanism  on  k,  is  indirect.  Its  analytical  expressions  are  practically  not  obtainable. 

Due  to  the  Robin  boundary  condition  type,  the  growth  rate  ifc,  is  embedded  in 
the  matrix  equation  that  determines  the  source  distribution  for  diffusion,  as  shown 
schematically  below: 

A(k,)a  =  b. 

It  requires  an  implicit  or  explicit  inversion  of  the  matrix  A  to  obtain  o.  Given  the 
largeness  of  A,  it  is  thus  impossible  to  maintain  analyticity  of  k,  in  the  solution  and 
perform  differentiation  on  it  later  to  get  the  Jacobian.  For  now,  assume  somehow  the 
Jacobian  can  be  obtained,  the  change  in  oxide  injected  into  bulk  due  to  a  sma.!!  change 
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in  k,  is  given  by 

AiI=kJDAk,.  (4.38) 

where  A;  is  an  appropriate  scale  factor. 

The  situation  with  oxide  flow  is  simpler.  The  relationship  between  the  injection  u 
and  P  is  linear  and  is  given  by 

P  =  AU  +  Po  (4.39) 

where  Po  is  contributions  from  nitride  or  free  surface  segments.  Hence  its  Jacobian  is 
given  by  Jp  =  A.  The  change  in  P  due  to  k,  is  therefore  given  by 

APn  =  JpJD^k,.  (4.40) 


Approximation  on  Diffusion  Jacobian 


As  mentioned  earlier,  the  diffusion  Jacobian  Jd  cannot  be  obtained  by  analytical 
means.  A  numerical  check  on  a  typical  oxidation  problem  has  show  that  it  can  be 
approximated  with  a  diagonal  matrix.  The  oxide  growth  rate  of  a  one-dimensional 
problem  is  given  by  Eq.  2.6: 


1  jk,C* 

N  k.k,do 

^  +  T  +  ~zr 

1  k,C* 

1  +  -^ 

Its  sensitivity  to  k,  is  defined  as 

dG  C*  £>* 
dk,~  N  (Jfe.do  +  Dy 


(4.41) 


When  the  oxide  is  thin,  the  growth  rate  is  reaction-rate  limited.  In  that  case,  ^ 
is  almost  C*/jV,  the  maximum  possible.  As  the  oxide  gets  thicker,  the  growth  rate 
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is  diffusion  limited  and  the  sensitivity  becomes  weak:  [C'fN)  •  [D^ Such  a 
prediction  does  not  generalize  to  a  two-dimensional  system.  An  increase  in  k,  of 
a  particular  boundary  location  does  not  necessitate  an  increased  supply  of  oxidant 
from  the  free  surface  to  enhance  its  growth  rate.  Rather,  it  can  absorb  oxidants 
that  originally  botmd  for  neighboring  points.  Thus  the  sensitivity  of  growth  rate 
on  k,  is  higher  than  the  one-dimensional  case  for  thick  oxide.  Through  numerical 
investigations,  it  has  been  found  to  be  close  to  C*/JV  in  most  cases.  More  surprisingly 
and  perplexing,  the  growth  rate  of  neighboring  points  are  not  ziffected  much.  Thus  a 
good  approximation  for  the  Jacobian  of  the  diffusion  problem  is  just  a  diagonal  matrix. 


Approximation  on  Jacobian  for  Oxide  Flow 

A  change  of  k,  in  one  segment  can  in  principle  affect  the  rest.  However,  due  to  ge¬ 
ometric  and  feedback  considerations,  strong  influence  is  limited  to  a  few  neighboring 
segments;  the  rest  experience  very  little  changes.  The  implication  on  the  Jacobian  is 
that  terms  far  away  from  the  diagonal  are  apprmcimately  0  (provided  the  k/s  are  prop¬ 
erly  ordered).  In  our  tridiagonal  approximation,  entries  that  are  not  on  the  diagonad 
or  the  two  off-diagonals  are  set  to  zero. 

The  oxide  flow  problem  is  not  subject  to  the  isolation  treatment  as  the  diffusion. 
Numerical  calculations  of  stress  due  to  the  increased  growth  rate  of  a  segment  yield  a 
— / +/—  pattern  distributed  over  a  few  segments.  Compressive  stress  is  experienced  by 
the  modifying  segment,  and  tensile  stress  by  its  two  adjacent  segments.  The  next  two 
further  out  also  see  tensile  stress,  but  with  much  less  magnitude.  In  consideration  of  the 
approximation  employed  previously  in  the  diffusion  system,  there  is  no  strong  reason 
to  model  stress  change  beyond  the  two  neighboring  segments.  Thus  the  approximated 
Jacobian  is  a  tridiagonal  matrix. 


The  remaining  problem  is  to  calculate  the  entries  of  the  matrix,  that  is  the  esti¬ 
mation  of  stress  due  to  a  unit  change  in  the  injection  rate  of  each  silicon  segment. 
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Suppose  among  the  N  boundary  segments,  M  of  them  are  for  the  silicon  interface. 
In  the  original  solution  scheme,  we  solve  for  the  source  S  in  terms  of  the  boundary 
condition  6: 

>15  =  6. 

The  method  is  extended  to  find  simultaneously  the  solutions  for  the  “impulse”  re¬ 
sponses  associated  with  the  M  silicon  segments.  The  single-col\imn  vectors,  5  and  6, 
are  expanded  to  (Af  +  l)  columns.  After  the  normal  displacement  components  from 
the  silicon  segments  are  removed,  the  modifed  b  is  put  in  the  zeroth  column  of  the 
new  multi-column  6.  Unit  vectors  corresponding  to  the  missing  normal  displacement 
component  are  placed  in  each  of  the  remaining  columns.  After  solving  for  5,  we  can 
reconstruct  the  desired  source  distribution 

s  =  s|oi+f;u„ii]S[.i. 

*=x 

5[1]  to  S[N]  are  unit-impulse  solutions  from  which  we  can  obtain  the  stress  contri¬ 
bution  due  to  each  silicon  segments,  that  b,  the  entries  for  the  Jacobian.  Although 
it  requires  significant  computation  time  to  obtain  5,  considerable  savings  are  later 
achieved  for  each  subsequent  iteration  because  the  matrix  equation  is  not  solved  again. 

Given  that  the  Jacobians  for  diffusion  and  oxide  fiow  are  a  diagonal  matrix  and 
a  tridiagonal  matrix  respectively,  the  combined  Jacobian  is  still  a  tridiagoned  matrix. 
Because  the  Jacobians  are  crudely  approximated,  the  parameter  Sk  is  utilized  to  make 
the  system  more  diagonally  prominent.  Typically  Sk  b  chosen  to  be  less  than  the 
largest  entry  in  J^Jd- 

Nonlinear  Diffusivity  and  Viscosity 

The  feedback  mechanism  for  the  bulk  parsimeters  diffusivity  and  vbcosity  b  more 
indirect  than  the  reaction-rate.  Consider  the  diffusion  problem.  Suppose  the  diffusivity 
of  a  region  b  perturbed  and  made  smaller  than  its  equilibrium  value.  Thb  reduces 
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the  oxidant  flxix  to  a  certain  section  of  the  silicon  interface.  The  oxide  growth  rate  at 
that  section  is  altered,  aoid  that  consequently  generates  a  new  stress  field  that  reaches 
the  perturbed  region  and  changes  its  diffusivity.  The  details  depend  on  how  far  the 
perturbed  region  is  from  the  silicon  interface.  The  situation  with  nonlinear  viscosity 
is  also  complicated.  Without  any  knowledge  of  their  Jacobizms,  these  two  parameters 
are  solved  using  the  simple  iterative  method. 


Chapter  5 
Simulations 


Simulation  results  for  various  geometries  and  conditions  are  presented  in  this  chapter. 


5,1  Viscoelastic  Stress  Relaxation 

We  will  consider  simple  cases  of  local  oxidation  to  demonstrate  how  stress  is  relieved 
through  viscoelastic  flow.  A  wide  range  of  stress  relaxations  times  axe  tested  to  see 
how  difi'erent  relaxation  speed  changes  the  stress  behaviors  eoid  distributions.  No  stress 
effects  are  included  in  this  study. 

Fig.'  5.1  shows  the  simulation  results  for  a  semi-recessed  LOCOS  structure  for 
different  stress  relaxation  times.  In  (a),  the  outline  of  the  structure  is  plotted  at  every 
time  step.  In  this  simulation  window  of  width  the  silicon  nitride  mask  extends 

from  I  =  0  to  z  =  0.96/im  and  is  assumed  to  be  totally  flexible.  The  pad  oxide 
thickness  is  200  A.  Oxidation  is  carried  out  at  900'’C  in  a  wet  ambient  for  4.6  hours  to 
obtain  a  flnal  fleld  oxide  thickness  of  5000  A.  Young’s  modulus  and  Poisson’s  ratio  are 
taken  to  be  8  x  10“  dynes-cm~*  and  0.194.  Shown  in  b,  c  and  d  is  stress  distribution 
corresponding  to  relaxation  times,  tj/G,  of  1  minute,  1  hour,  and  10  hours  respectively. 
The  normal  component  of  the  surface  traction  at  the  oxide-silicon  interface  is  plotted 
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Figure  5.1:  The  nitride  mask  is  assumed  to  be  totally  flexible,  (a)  Boundary  outline. 
Stress  behavior  for  stress  relaxation  time  of  (b)  1  min;  (c)  1  hour;  and  (d)  10  hours. 
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at  every  time  step,  like  the  outline  of  the  oxide.  In  all  three  stress  plots,  there  are 
two  peaks  in  the  compressive  stress  region  (the  negative  value  fange).  The  early  peak 
occurs  at  the  edge  of  the  nitride  mask.  This  peak  is  due  to  the  highly  nonuniform 
oxidation  rate  in  that  region.  As  time  progresses,  the  peak  shifts  to  the  left,  further 
under  the  nitride  mask,  and  gives  rise  to  a  late  peak.  As  one  would  expect,  stress 
reduces  as  viscosity  decreases.  Maximum  stress  attained  for  the  3  relaxation  times  is 
8.1  X  10®,  1.5  X  10^°,  and  2  x  10^°  dynes-cm~*  respectively.  The  oxide  shape  is  not 
affected  much  by  the  relaxation  times,  therefore  only  one  is  shown. 

We  now  repeat  the  same  simulations  with  the  silicon  nitride  mask  modeled  as  an 
elastic  material.  The  thickness  of  the  nitride  layer  is  O.l^m.  The  final  shapes  of  the 
oxide  and  nitride  layer  are  shown  in  Figs.  5.2a,  b,  and  c.  As  the  oxide  becomes  less 
viscous  and  flows  more  readily,  the  nitride  layer  bends  less.  The  corresponding  values 
for  peak  stress  are  3.8  x  10^°,  1.2  x  10^^,  and  1.6  x  10^^  dynes-cm~*  respectively.  The 
last  two  vaJues  are  unrealistically  high  as  they  are  of  the  same  order  of  magnitude  as 
Young’s  modulus  of  oxide  or  nitride.  We  have  found  that  in  general  we  need  these 
high  stress  values  in  order  to  deflect  the  nitride  mask  to  a  degree  comparable  to  what 
is  found  in  experiments.  To  keep  stress  down  to  a  realistic  level,  the  nitride  mask 
must  deform  plastically  or  viscoelastically.  Unfortunately  we  don’t  have  any  models 
for  those  behaviors  nor  any  data  to  determine  a  model.  In  all  likelihood,  the  stress 
value  is  inaccurate  if  it  is  a  few  percent  of  Young’s  modulus. 

Viscoelasticity  is  a  relative  concept.  If  the  time  period  involved  is  much  shorter 
than  the  time  constant  of  stress  relaxation,  a  viscoelastic  material  deforms  elastically. 
If  the  converse  is  true,  the  material  flows  viscously.  The  viscoelastic  behavior  is  most 
profound  when  the  oxidation  time  and  relaxation  time  are  comparable.  Shown  in 
Fig.  3a  is  a  plot  of  maximum  stress  verstis  the  relaxation  time.  The  simulation  condi¬ 
tions  are  the  ones  used  earlier,  without  modeling  the  elastic  bending  of  nitride  mask. 
For  long  relaxation  times  (r  >  10  hours),  the  stress  curve  is  constant  as  the  mate¬ 
rial  is  essentially  elastic.  For  short  relaxation  times  (r  <  1),  the  curve  bends  down 
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Figure  5.2:  The  nitride  is  assumed  to  be  elastic.  Final  oxide  profile  is  shown  for  oxide 
relaxation  time  of:  (a)  1  min;  (b)  1  hour;  and  (c)  10  hours. 
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with  a  constant  slope.  The  drop  is  a  result  of  stress  being  inversely  proportional  to 
the  relaxation  time.  For  all  the  data  points,  maximum  stress  is  observed  at  around 
35  minutes  into  oxidation;  that  value  roughly  corresponds  to  the  transition  region  in 
Fig.  5.3a.  An  alternate  way  of  assessing  the  significamce  of  viscoelastic  flow  is  to  check 
the  sensitivity  of  stress  to  a  change  in  viscosity.  If  the  oxide  is  elastic,  stress  will  not 
be  affected.  On  the  other  hand,  if  it  is  viscous,  one  would  see  the  same  change  in 
stress.  Shown  in  Fig.  5.3b  is  the  stress  change  due  to  a  10%  increase  in  viscosity. 
There  is  little  change  in  stress  value  for  large  relaxation  times  and  an  almost  full  10% 
change  for  short  relaxation  times.  This  essentially  corresponds  to  the  gradient  of  the 
cTirve  in  Fig.  3a.  As  we  can  see,  the  stress  value  levels  off  after  the  viscosity  reaches 
a  certain  threshold.  If  a  viscous  flow  model  was  used  instead,  stress  would  be  strictly 
proportional  to  the  viscosity. 


5.2  Nonhomogenous  Bulk  Parameters 


Here  we  demonstrate  the  capability  of  the  initial  stress/built-in  field  formulation  for 
modeling  stress  effects  on  viscosity  and  diffusivity.  Concave  md  convex  silicon  comers 
are  oxidized  at  950^Cin  a  wet  oxidation;  approximately  1100  Aof  oxide  is  grown. 

Fig.  5.4  shows  the  oxidation  of  a  convex  silicon  comer  and  the  domain  partition 
used  for  computing  pseudo  domain  sources.  When  the  stress  effect  on  viscosity  is  not 
modeled,  the  peak  stress  goes  as  high  as  2.7  x  10^^  dynes/cm“*.  With  a  stress  model 
included,  the  peak  stress  drops  down  to  2.2  X  10®  dynes/cm“®. 

The  retardation  effect  of  stress  on  diffusivity  is  illustrated  in  Fig.  5.5.  Oxidation  is 
done  on  an  concave  silicon  comer.  A  large  compressive  hydrostatic  pressure  is  present 
at  the  comer.  In  (a) ,  no  stress  effect  on  diffusivity  is  considered  -  no  growth  retardation 
at  the  comer  is  observed.  In  (b),  stress  effect  is  modeled.  The  oxide  thinning  observed 
at  the  comer  is  due  to  reduced  diffusivity  of  oxidant. 


Percentage  Change  In  Stress  (%)  Maximum  Stress  (dynes/cm*2) 
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Maximum  Stress  versus  Relaxation  time 


Stress  Change  versus  Relaxation  time 


Relaxation  Time  (hrs) 


(b) 


Figure  5.3:  (a)  Maximum  stress  versus  relaxation  times,  (b)  Stress  change  due  to  a 
10%  increase  in  viscosity  at  different  relaxation  times. 
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Fig:ure  5.4:  Modeling  Stress  Effect  on  Viscosity,  (a]  shows  the  outlines  of  oxide  on  a 
convex  silicon  comer,  (b)  shows  the  domain  partition  for  modeling  nonlinear  viscosity. 
Cells  are  placed  at  the  comer  where  high  stress  is  expected. 
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5.3  Modeling  Kao’s  Experiments 


In.  this  section,  we  demonstrate  a  simplified  viscoelastic  fiow  model  for  fitting  Kao’s 
experimental  data.  This  model  ignores  the  crystallographic  dependency  of  the  re¬ 
action  rate  k„  it  solves  radial  symmetric  problems  whose  solutions  are  pseudo  one¬ 
dimensional  in  nature.  The  oxide  is  assumed  to  be  viscoelastic.  An  effective  diSusivity 
and  viscosity  is  assumed  for  the  domain.  Shown  below  are  the  stress  models: 


k, 

D 


<TYVn2lkT 

sinh(<rrV’,2/ikr) 


where  ay  —  {Cy  is  the  maximum  distortion  energy.)  No  maximum  limit  is 

placed  on  any  parameters.  With  the  following  fitting  values, 
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reasonable  simulation  results  have  been  obtained.  They  are  compared  with  Sutardja’s 
in  Fig.  5.6.  As  it  can  be  seen,  they  all  agree  very  well.  However,  the  fittings  values  are 
not  close.  For  instance,  our  Vj,  and  Vd  are  45  and  50  A®  respectively,  as  compared 
to  12.5  A’  and  75  A^  in  Sutardja’s  model.  Besides  the  differences  in  formulation  and 
numerical  approximations,  the  disagreement  may  also  be  attributed  to  the  scaling 
problem  of  activation  volumes,  which  b  explained  in  the  next  section. 
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5.3.1  Arbitrary  Scaling  of  Activation  Volumes 


We  have  more  unJcnown  parameters  than  data  points  can  fit  for  thermal  oxidation 
processes.  Part  of  the  problem  is  that  we  are  imable  to  measure  local  stress  m  situ 
or  other  activities.  Experimental  measurements  are  usually  given  in  terms  of  oxide 
thickness  -  a  collective  parameter  due  to  many  stress  effects.  It  is  difficult  to  come  up 
with  a  imique  set  of  activation  volumes  for  the  stress  effects. 


A  problem  inherit  in  incompressible  viscous  flow  methods  is  that  the  activation 
volumes  for  stress  parameters  are  scalable  -  their  values  can  be  changed  without 
altering  the  magnitudes  of  retardation  effects.  The  root  of  the  problem  is  that  stress 
is  proportioned  to  the  viscosity  and  that  all  the  stress  equations  are  of  the  form 


which  is  restructmed  to  show  the  dependency  on  t;: 


we  realize  that  one  important  consideration  in  this  expression  is  the  product  consisting 
of  viscosity  and  activation  volume,  riV.  As  long  as  this  product  is  kept  constant, 
the  value  of  /  will  be  the  same,  and  hence  the  stress  effects  will  be  the  same  too. 
Hence,  doubling  all  activation  volumes  and  halving  the  viscosity  will  produce  the  same 
reduction  in  oxide  thickness,  although  the  stress  values  are  reduced  by  half.  We  are 
free  to  choose  any  viscosity /activation  volumes  combination  to  get  any  stress  values 
we  want.  Kao  eliminates  this  degree  of  freedom  by  fixing  Vi,,  the  activation  volume  for 
to  be  25A,  which  is  the  net  volume  increase  of  an  oxidized  Si  atom.  In  Sutardja’s 
model,  Vi,  is  12.5A.  Probably  it  has  to  do  with  the  silicon  monoxide  phrase  that  an 
Si  atom  has  to  go  through  in  the  oxidation  event.  In  any  case,  unless  we  have  better 
measurement  of  oxide  viscosity  or  stress,  the  values  for  activation  volumes  will  still  be 
debatable. 


We  have  also  attempted  to  model  Kao’s  experiments  in  true  two  dimensions,  taking 
in  accoimt  the  orientation  dependency  of  the  reaction  rate.  Unfortunately,  the  results 


Figure  5.6:  Simulations  of  Oxide  Growth  on  Cylindrical  Structures,  (a)  Simplified 
viscoelastic  flow  model,  (b)  Incompressible  viscous  flow  model. 
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do  not  matched  very  well  with  the  experimental  data.  Possibly,  part  of  the  problem  has 
to  do  with  the  fact  that  stress  history  is  not  propagated  from  one  timestep  to  another 
when  operating  under  nonlineax  viscosity  condition,  due  to  progreiming  constraints. 


Chapter  6 


Conclusions 


Thermal  oxidation  of  silicon  involves  two  parallel  processes,  namely  the  motion  of 
oxide  caused  by  volume  expansion  and  the  diffusion  of  oxidmts  such  as  H2O  or  O2. 
An  “indirect”  boundary  element  approach  for  modeling  oxidation  has  been  introduced. 
The  diffusion  problem  is  modeled  as  a  standard  Laplace  problem  using  the  scalar 
potential  formulation;  the  motion  of  oxide  is  solved  with  a  vector  potential  formulation. 

Three  different  motion  models,  namely,  viscous  flow,  elastic  deformation,  and  vis¬ 
coelastic  flow,  are  considered.  The  BEM  formulation  for  these  models  are  essentially 
identical,  with  the  only  difference  in  the  kernel  functions  (or  Green’s  functions).  A 
Laplace  transform  technique  is  used  to  trsmsform  Kelvin’s  solution,  which  is  a  funda¬ 
mental  solution  for  elasticity,  into  a  viscoelastic  kernel  function.  A  constant-velocity 
loading  function  is  chosen  to  ensure  a  wide  range  of  operation  conditions.  Essentially, 
it  can  model  elastic  deformation  as  well  as  incompressible  viscous  flow  problems. 

Stress  generated  during  oxidation  steps  affects  many  variables  -  some  boundary 
parameters,  some  bulk  parameters.  Diffusivity  of  oxidants  and  viscosity  of  oxide  are 
the  bulk  parameters  rendered  nonhomogeneous  and  nonlinear  by  stress.  For  them, 
domain  calculations  are  needed.  A  imifled  method  for  dealing  with  nonhomogeneity 
in  both  diffusion  and  viscoelastic  flow  has  been  developed.  This  initial  stress/built- 
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in  field  scheme  utilizes  pseudo  domain  sources  to  generate  correction  terms  to  the 
solutions. 


6.1  Future  BEM  in  Thermal  Oxidation 


For  what  it  is  designed  for,  namely,  linear  systems,  the  BEM  performs  very  well. 
In  fact,  it  has  been  found  to  be  superior  to  the  finite  element  method  in  fracture 
mechanics.  In  applications  where  stress  is  highly  concentrated  in  a  small  region,  the 
BEM  is  able  to  generate  good  results  without  much  refinement  on  the  boundary  seg¬ 
ments.  In  contrast,  the  finite  element  method  requires  fine-tuning  of  the  mesh  to  get 
accurate  results.  The  superior  performance  of  the  BEM  comes  from  Green’s  function 
that  is  singular  -  it  is  difficxilt  to  deal  with,  but  it  can  model  large  change  in  stress 
field  accurately. 

The  situation  with  nonlinear  problems  is  less  clearly  defined.  In  dealing  with  non- 
homogeneous  materials,  the  BEM  requires  domain  calculations  with  either  interior  cells 
or  subregions.  In  places  where  strong  variations  in  the  bulk  p2Lrameters  are  observed, 
large  number  of  interior  cells  or  subregions  may  be  required.  Granted,  no  cells  or 
subregions  are  needed  in  regions  where  stress  is  insignificant.  However,  for  high  stress 
areas,  the  cell  generation  requirements  may  be  simile  to  those  of  the  FEM. 

Given  the  fact  that  the  matrix  setup  operation  is  O(n^)  and  the  solution  is  O(n^), 
where  n  is  total  number  of  boundary  elements  and  interior  cells  (or  interfaciaJ  sources), 
a  severe  computational  penalty  is  imposed  on  the  BEM  when  there  is  a  sizable  popu¬ 
lation  of  cells. 

To  determine  the  viability  of  BEM  in  thermal  oxidation  modeling,  more  inves¬ 
tigative  works  are  needed  on  iterative  methods  for  nonlinear  solutions.  Currently,  the 
domain  pseudo  sources  and  boundary  sources  are  solved  simultaneously.  In  other  BEM 
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applications,  the  domain  perturbation  is  not  incorporated  in  the  system  matrix  but 
calculated  externally  in  an  iterative  manner.  Such  schemes  help  to  reduce  the  matrix 
solution  time  significantly. 


6.2  Proposed  Experimental  Studies 


In  order  to  function  effectively,  process  simulators  need  accurate  models  amd  pauameter 
values.  To  determine  the  effects  of  stress  on  various  parameters  in  am  oxidation  model, 
stress  must  be  introduced  to  the  test  samples  in  a  controlled  manner  so  that  the  results 
can  be  quantified.  So  far,  Kao’s  experimental  technique  of  using  cylindrical  structures 
represents  one  of  the  best  approaches. 

More  experimental  studies  and  measurements  can  be  done  with  cylindrical  struc¬ 
tures.  Consider  the  crystal  orientation  on  the  reaction  rate.  Only  the  thinning  effects 
on  the  [110]  surfaces  have  been  presented.  The  [lOO]  surfaces  should  be  able  to  provide 
another  set  of  data  for  simulation  verifications.  Before  that  is  done,  we  need  to  know 
how  the  reaction  rate  changes  as  the  crystal  orientation  varies  from  [110]  to  [lOOj.  This 
is  important  because  stress  distribution  is  sensitive  to  the  orientation  dependency  of 
the  reaction  rate.  Different  oxide  thicknesses  should  also  be  tried  to  extract  more 
information.  Consider  the  stress  behavior.  Pressure  is  primarily  a  function  of  the 
free  surface  curvature;  however,  the  surface  traction  on  the  silicon  interface  depends 
strongly  on  the  thickness  of  the  oxide.  Thus,  growing  different  oxide  thicknesses  may 
help  to  differentiate  the  contributions  of  these  two  mechanisms  to  the  growth  retar¬ 
dation  phenomenon.  The  activation  volumes  in  stress  functions  are  scalable  for  an 
oxidation  model  based  on  incompressible  viscous  flow  of  oxide.  Unless  stress  or  vis¬ 
cosity  is  measured  independently,  there  will  be  uncertainty  in  the  range  of  values  for 
activation  volumes. 


Appendix  A 

Miscellaneous  Formulations 


The  formulae  for  principal  stresses,  stress  transformation,  and  analytical  line  integra¬ 
tions  axe  given  in  this  appendix. 


A.l  Principal  Stresses 


The  values  for  the  principal  stresses  are  given  by  the  roots  of  this  determinant: 

<Tli  —  <7  <7i2  Ciz 

O^ji  C22  —  <7  <^2Z 

<731  <732  <733  —  (7 

In  plane  strain,  <713  =  <733  =  <732  =  <731  =  0.  Thus,  the  determinzuit  is  reduced  to  the 
following  equation: 

-  <733)  [{<^11  -  cr){a22  -ff)  -  <7i2<72i]  =  0. 

The  roots  are  given  by: 
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They  are  ordered  so  that  cj  >  an  >  <7///. 
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A. 2  Stress  and  Strain  Transformation 


It  is  often  needed  to  transform  stress  and  strain  field  from  the  global  coordinate  frame 
to  a  local  frame.  The  expression  for  two-dimensional  stress  transformation  is 
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where  0  is  the  angle  measured  counterclockwise  from  the  x  axis  of  the  old  frame  to 
the  x'  axis  of  the  new  frame.  For  strain,  the  transformation  rule  is  slightly  different: 
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For  comparion,  the  rotation  expression  for  a  vector  is  show  here: 


smB  1  f  /i  1 
cosB  J  \  A  J 


(A.4) 


With  reference  to  Eq.  A.2,  one  notes  that,  for  the  uniaxial  stress  problem,  a  maximum 
<^12  occurs  when  9  =  45”,  as  mentioned  in  Chapter  2. 


A. 3  Special  Analytical  Line  Integrals 


In  determining  the  infiuence  of  a  source  distribution  on  the  collocation  point  of  the 
same  segment,  analytical  technique  must  be  employed  as  we  are  integrated  over  a 
singularity.  The  results  are  simple  because  the  collocation  point  is  located  at  the 
center  and  the  source  distribution  is  assumed  to  be  constant.  The  integration  path  is 
chosen  such  that  the  source  does  not  lie  on  the  boundary  but  just  outside  the  domain 
(0"^).  Shown  below  are  the  results  for  the  seven  basis  functions 
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where  I  is  the  length  of  the  segment,  and  0  the  tilt  angle  of  the  segment  with  respect 
to  the  xi  axis,  measured  counter-clockwise. 


A. 4  General  Analytical  Line  Integrals 


In  this  section,  we  consider  the  integral  results  for  the  case  when  the  location  of  the 
collocation  point  is  arbitrary  with  respect  to  the  source  segment.  For  simplicity,  the 
solutions  are  given  in  terms  of  a  local  reference  frame  in  which  the  source  segment 
(t.e.  the  integration  path)  appears  to  orientate  vertically  (in  the  y  direction),  going 
from  (ao,6o)  to  (^k),6i)* 
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